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PART  I 


INTRODUCTION 

The  form  of  the  differential  equations  used  to  describe  a  physical 
process  can  be  specified  from  basic  conservation  principles.  Very  often  the 
parameters  associated  with  the  differential  equations  are  unknown.  They  are 
determined  from  a  comparison  of  experimental  measurements  of  the  process  and 
the  solutions  of  the  differential  equations  that  describe  the  process. 
Parameter  estimation  is  an  integral  part  of  the  analysis  of  experimental 
data  in  terms  of  a  model  of  known  form  with  unknown  coefficients. 

There  has  been  very  little  work  done  in  the  area  of  parameter  esti- 
mation. The  parameters  or  coefficients  cannot  be  measured  directly.  The 
measurable  variables  are  generally  the  dependent  variables  of  the  differ- 
ential equations.  Therefore  it  is  not  simple  to  identify  the  parameters. 

The  estimation  of  parameters  in  ordinary  differential  equations  has 
received  considerable  interest  in  recent  technical  literature. 

Rosenbrock  and  Storey  [19]  have  presented  a  review  of  the  techniques 
of  generating  and  analyzing  parameter  estimates. 

Bellman,  Kagiwada  and  Sridhar  [20]  have  shown  how  techniques  from  non- 
linear filtering  and  estimation  theory  can  be  applied  to  the  estimation  of 
parameters  in  ordinary  differential  equations. 

Lee  [11]  has  presented  the  problem  of  parameter  estimation  using  the 
quasi! inearization  method. 

The  objective  of  the  present  work  is  to  estimate  parameters  in  the  dif- 
ferential equations  resulting  from  tubular  flow  chemical  reactors  with 


axial  diffusion.  The  physical  process  is  assumed  to  be  represented  by  a 
nonlinear  ordinary  differential  equation.  Linearization  is  carried  through 
by  considering  the  first  two  terms  in  the  Taylor's  series  expansion  of  the 
original  nonlinear  equation.  This  technique  is  a  generalized  Newton- 
Raphson  formula  for  functional  equations.  It  is  more  commonly  known  as 
the  quasi! inearization  method.  The  main  advantage  of  this  technique  is  that 
the  procedure  converges  quadratically  to  the  solution  of  the  original 
equation,  if  it  should  converge. 

Integration  of  the  linearized  equation  is  carried  through  by  the  fourth 
order  Runge-Kutta-Gill  method.  An  algorithm  is  devised  based  on  the  least 
squares  method.  Initial  guesses  are  made  for  the  parameters  and  also  for 
the  function  value  of  concentration  as  a  function  of  axial  length.  Initial 
values  for  the  numerical  integration  are  so  chosen  that  the  integration  con- 
stants are  identically  equal  to  the  unknown  parameters.  Parameters  thus 
obtained  are  used  recursively  to  obtain  improved  results. 

Errors  in  the  experimental  data  are  considered  small  in  magnitude. 
The  mathematical  model  is  devised  to  fit  the  data. 


CHAPTER  I 
QUASILINEARIZATION  TECHNIQUE 

Solutions  to  initial  value  problems  are  well  developed  theoretically. 
They  are  also  easily  adaptable  for  solving  on  high  speed  digital  and  analog 
computers.  Nonlinear  differential  equations  having  two  or  multipoint 
boundary  values  are  common  in  engineering  and  physical  sciences.  There  is 
no  general  proof  for  the  existence  and  uniqueness  of  the  solutions  of  such 
problems.  Numerical  difficulties  in  their  solution  are  caused  because  not 
all  the  conditions  are  given  at  one  point.  A  trial  and  error  procedure  is 
generally  used  to  obtain  the  missing  initial  condition.  This  procedure  has 
a  relatively  slow  rate  of  convergence. 

Quasi! inearization  is  a  useful  technique  for  obtaining  numerical  solu- 
tions for  these  type  of  problems.  In  this  method  the  nonlinear  differential 
equation  is  solved  recursively  by  a  series  of  linear  differential  equations. 
The  quasilinearization  technique  converges  quadratically  to  the  exact 
solution,  if  at  all  it  should  converge.  Quadratic,  convergences  implies 
that  the  error  in  every  succeeding  iteration  tends  to  be  proportional 
to  the  square  of  the  error  in  its  immediately  preceeding  iteration. 

Linearization  of  the  original  differential  equation  is  carried  through 
by  considering  the  first  two  terms  in  the  Taylor's  series  expansion.  This 
is  a  generalized  Newton-Raphson  technique  for  functional  equations.  The 
Newton-Raphson  technique  is  associated  with  two  important  properties. 
These  are  monotone  convergence  and  quadratic  convergence.  The  nature  of 
the  monotone  convergence  property  depends  on  the  function  itself.  The 


monotone  convergence  property  exists  for  the  Newton-Raphson  formula  only 
if  the  function  is  a  monotone  decreasing  or  monotone  increasing  function. 
The  function  should  be  strictly  convex  or  strictly  concave.     In  general   the 
Newton-Raphson  formula  always  has  the  quadratic  convergence  property  pro- 
vided that  it  converges  [11]. 

The  Newton-Raphson  equation  is  always  linear  even  if  the  original 
function  is  nonlinear.     Linear  boundary  value  problems  can  be  solved 
fairly  routinely  on  modern  computers. 


CHAPTER  II 

LINEARIZATION  OF  THE  NONLINEAR  ORDINARY  DIFFERENTIAL  EQUATION 

The  mathematical  formulation  for  estimating  the  parameters  in  the 
differential  equation  for  homogeneous  tubular  flow  chemical  reactor  with 
axial  mixing  is  presented  in  this  chapter.  The  physical  process  is  assumed 
to  be  represented  by  a  nonlinear  second-order  ordinary  differential  equation. 

1  d2x   dx_Rx2=Q  (1) 


Pdt2   dt 


where, 


Lv 
P  is  the  dimensionless  Peclet  group.  Its  magnitude  is  -jj-  , 

kL 
R  is  the  reaction  rate  group  —  , 

t  is  the  dimensionless  reactor  length  which  varies  between  0  and  1. 
It  is  obtained  by  dividing  the  actual  position  along  the  axial 
direction  of  the  reactor  by  the  total  reactor  length  L, 

x  is  the  concentration  of  the  reactant, 

v  is  the  flow  velocity  of  the  reaction  mixture.  It  is  assumed 
constant  throughout  the  reactors, 

D  is  the  mean  mass  axial  dispersion  coefficient.  It  is  assumed  con- 
stant, 

k  is  the  specific  chemical  reaction  rate.  It  is  also  assumed  con- 
stant. 


Boundary  Conditions 

xe  -  x(0)  -  14|  |t=Q  =  C)     say  at  t  =  0  (2) 


dx 


=0  at  t  =  t,.     ,  (3) 

t=tf  final  ' 


where , 


x„  is  the  concentration  of  the  reactant  before  it  enters  the  reactor, 
e 

It  is  a  known  quantity,  x(0)  is  the  concentration  just  after  entering 
the  reactor.     There  is  a  discontinuity  of  the  concentration  x  of  the  re- 
actant at  the  entrance  to  the  tubular  reactor. 

Equation  (1)  is  a  second  order  nonlinear  differential   equation  of 
the  boundary-value  type  with  P  and  R  as  parameters.     These  parameters 
are  unknown  values  that  are  to  be  determined. 

The  preliminary  step  involved  for  estimating  the  parameters  is  to  lin- 
earize the  original  nonlinear  equation  (1)  by  the  generalized  Newton- 
Raphson  method. 

£f=p£+PRx2  (5) 

Equation  (5)  can  be  rewritten  as, 

x'  =£=Y  (6a) 

y'  =  %=  Py  +  PRx2  (6b) 


The  parameters  P  and  R  are  considered  as  dependent  variables  parallel 
to  x  and  as  functions  of  the  independent  variable  t.  The  parameters 
satisfy  the  following  equations, 


and 


f-0.  (60 


f-o.  («) 


The  generalized  Newton-Raphson  formula  is  applied  to  equation  (6a)  as 

follows. 

8x' 
Xn+1  =  xn  +  ^7  {Vl  "  V  (7a> 

j_  L  —J- 

The  subscripts  n  and  n+1   denote  the  n      and  (n+1)       iteration  respectively. 
From  equation  (6a) , 
x'  =  y. 
Thus, 

xn+l   =^n+ly7{Vl   "  V  (7b) 

xn+1  -  yn  +  l  •  (yn+1  -  yn}  (7c) 

xn+l   "  Vl  (8a) 

Similarly  equation  (6b)   is  written  as, 

y'  =  Py  +  PRx2  (7d) 

Thus, 

y1  =  f(x,y,P,R)  (7e) 


^1   =  yn  +  ^  (Vl   "  V  +  5*  ^   '  ^ 
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sy' 


9y' 


+  l^Pn+1    ■  V+^(Rn+l   -Rn) 


(7f) 


Thus, 


2s         3 


*rt  =  <  Vn  +  PnRnV  +  aF  <Vn  +  PnRnxn><Vl  "  xn> 


+  W  (Vn  +  Wn^n-H  "  V  +  W~  <  Vn  +  PnRnxn^Vl  '  Pn> 


+  4(Pnyn  +  P"RnX")(Rn+l   "  Rn> 


(7g) 


+  (Pn-   l-0)(yn+1   -yn)Ml'y/l-Vn)(Vl   -Pn) 


+  (0  +  Pn  .   1    .  <)(R 


"  Rn) 


(7h) 


*ifl  =  (pnyn  +  PnVn>  +  (2PnRnxn>(xn+l  "   xn> 

+  Pn(yn+1  -  yn)  +  (yn  -  V>n+i  -  Pn) 


+  PnVRn+l   "  V 


(71) 
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Rearranging  the  terms  in  equation  (7i) 


*rt  -  PnVl  +  2PnRn*nVl  +  VnRn+l  +  <*n  +  Vn'Vl 


(3Pn>Vn  +  W  <8b' 


Similarly  equations  (6c)  and  (6d)  are  written  as, 


R^+1  =  0,  (8c) 


The  boundary  conditions  (2)  and  (3)  are  written  below  as 

xe  =  x(0)  -  p-y(O)  =  c  at  t  =  0  (9a) 

y(tf)  =0  at  t  =  tf       •  (9b) 

x(t-,)  =  b,  at  t  =  tj  (9c) 

x(t2)  =  b2  at  t  =  t2  (9d) 


where  b-,  and  b2  are  the  true  values  of  x  at  t,  and  t2. 

Applying  the  generalized  Newton  Raphson  method  to  equations (9), 


xn+1(0)  =  c  (10a) 

yn+1(tf)  =  o  (iob) 
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Vi<V  =  bi  '  (10c> 

Vl«V  -  b2  (10d) 

Equations  (8)  together  with  the  boundary  condititions  (10)  are  solved 
by  the  use  of  the  principle  of  superposition.  Since  one  initial  condition 
equation  (2)  is  given,  three  sets  of  homogeneous  solutions  are  needed  for 
solving  the  equations  (8).  Initial  values  for  solving  equations  (8)  are 
so  chosen  that  they  satisfy  the  initial  boundary  condition  and  also  set  the 
integration  constants  equal  to  the  unknown  parameters. 

The  particular  solution  is  obtained  with  the  following  initial  values. 

xpin+1(0)=c  (11a) 

W°>  =  °  (llb) 

RP;„+1(0)=0  (lie) 

wo)-°  (lld) 

where,  the  subscript  p  denotes  the  particular  value. 

The  homogeneous  solutions  are  obtained  by  integrating  the  following 
equations. 

K*\  ■  Vi  (12a) 


y;+1  ■  p„yn+1  ♦  2PnV„Vi  +  VnVi  +  <y, +  V>n+i  "2b» 
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r;+1  =  o  (i2c) 

P.+1  .  o  (12d) 

The  three  sets  of  initial  values  of  obtaining  the  homogeneous  solution 
are  given  below.     Subscripts  h-, ,   h„,  h-  denote  the  three  homogeneous  values. 

xhUnt,(0)  =  0  (13a) 

yhl,n+1(o)-i  03b) 

Rhl,n+1(o)  ■  0  (130 

Phl,n+1(0)-0  03d) 

xh2jn+1(0)  -  0  (14a) 

yh2,n+,(o>  =  0  (14b) 

Rh2in+1(0)  •  1  (140 

Ph2jn+1(0)-0  (140 


and, 


xh3in+1(0)  -  0  (15a) 

yh3>n+1(o)  =  o  (15b) 
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Rh3,n+1(0)  -  0  (lie) 

Wl(0>  =  1  (15d) 

The  general   solutions  of  the  system  of  equations  (8)  are  obtained  by  the 
principle  of  superposition. 


W*>-'Wl<t>+-J1   aJ>n+l   *hi.«*W  (lfa) 


3 

i 


VlW'Vn+l**'*!   aj,n+l  *hj,n+l(t)  (16b) 


■W*>*  Wt,  +  j2,   a0,n+l  «M.i»l(t)  (16C) 

W*>"  Wt}*i   aj,n+l   Wl(t>  <16d> 

J  ' 

It  will   be  seen  from  equations  (16c)  and  (16d)  that, 

Rn+l(°)=a2sn+l  (17) 

Vl(°)  =  a3,n+1 

It  was  initially  assumed  that  both  Rn+1 (t)  and  Pn+1(t)  are  constant 
functions.  Thus  equations  (17)  and  (18)  are  true  for  all  values  of  t. 
Thus, 

W*>  ■■2.1*1  (17a) 


(18) 
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Pn+l(t>  =  a3,n+l 
The  integration  constant  a,  n+1  can  be  expressed  as  a  function  of 


(18a) 


a        n  and  a0      ,.     The  following  equation  is  obtained  from  equation  (16b) 
2.,  n+1  3,n+l 


at  t  =  tf 


a 


frp.n+^V  *  a2,n+l  yh2,n+l(tf}  +  a3.n+l  yh3,n+l(tf)}     (19) 
l.n+1  *hl,n+llV 
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CHAPTER  III 
NUMERICAL  ANALYSIS 

The  fourth  order  Runge-Kutta  Gill  method  has  been  used  for  carrying 
through  the  numerical   integration.     Integration  over  the  interval   (0,1) 
is  carried  through  with  a  step  size  of  0.01. 

Assuming  that  a  solution  exists,  the  linear  differential   equation  is 
solved  in  two  steps.     First  one  set  of  particular  and  three  sets  of  homo- 
geneous solutions  are  obtained  numerically.     Initial   values  for  obtaining 
these  have  been  discussed  in  the  preceeding  chapter.     Then  the  integration 
constants  are  obtained  by  the  least  square  method. 

The  error  between  the  general  solution  and  the  experimental  values 
is  minimized  as  follows 

m1 

IU-    J    [*n+1(ts)-bs]2  (1) 


ml 


(t.)  + 


Qn+1  =     I    [xp,n+l(ts)  +  al,n+l  xhl,n+l(ts}  +  a2,n+l  xh2,n+l(ts 
s=l       r' 

a3,n+l   xh3,n+l(ts)  "  bs]  {2) 

From  equation  (19)  of  the  preceeding  chapter, 


a 


^pWV  +  a2,n+l  ^WM  +  a3,n+l  j^jjjVj     (3) 
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Introducing  equation  (3)  in  equation  (2), 


"'  *hl,n+1<V 


Qn+1  =  SI,   ^I'V  "  y^TV  <WV  +  a2.n+l  *h2.n+l<V 

+  a3,n+l  "h3.n+l(tf»  +  a2,n+l  xh2,n+l(ts>  +  a3,n+l  "hWV  "  bs]2 

(4) 

xhl   n+1     s 
Rearranging  the  terms  and  writing  A_,n   =  - — l  ■■■-(+  \ 

n   '       yh1,n+rV 
ml 

Qn+1  =  s|i  Cxpin+1(ts)  -  An+1  yp>n+1(tf)  +  a2jn+1(xh2>n+1(ts) 

"  Vl  yh2,n+l^f^  +  ^n+lKsWV  "  An+lyh3,n+l  (tf}>  "  bs]2 

(5) 

Since  all  the  particular  and  homogeneous  solutions  are  known  quantities 

equation  (5)  can  be  written  as 

m, 

Qn+1  =  J1   ^l,n+l(V  +  a2,n+l   '  ^2,n+l  (ts>  +  a3,n+l    '  t»3,n+l(ts)  "  bs]2 

(6) 


Minimizing  equation  (6)  wrt  a~  n+1  and  a3  n+1  the  following  two 
equations  are  obtained. 

ml 
^  iWV&l.irntV  +  a2,n+l    '  "2WV  +  a3,n+l    '  q3,n+l(ts^ 

-  bl  =  0 


(7a) 
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mi 
J,   ^WV^lWV  +a2,n+l    '   c'2,n+l(ts)  +  a3,n+l    *   q3,n+l(V 

-  bs]  =  0       (7b) 

The  two  simultaneous  algebraic  equations  (7a)  and  (7b)  are  to  be  solved 
to  obtain  the  integration  constants,  which  are  equal  to  the  unknown  para- 
meters. The  parameters  are  used  recursively  to  obtain  improved  results. 

The  flow  chart  figure  (2),  and  the  computer  program  for  solving  the 
problem  appear  in  the  Appendix. 
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CHAPTER  IV 
CONCLUSION 

The  quasi! inearization  method  combined  with  the  superposition  method  is 
unstable  under  certain  conditions  £11].  For  example,  during  the  process  of 
iteration  one  or  more  values  of  the  particular  and  homogeneous  solutions  can 
become  unreasonable.  The  solutions  will  not  converge  to  the  exact  values 
even  if  the  exact  values  of  the  parameters  are  used  as  the  initial  approx- 
imations. 

Extremely  large  or  small  values  of  the  parameters  are  obtained  at  the 
end  of  the  first  few  iterations.  Since  the  original  nonlinear  equations  are 
very  sensitive  to  these  values,  extremely  large  positive  or  negative  values 
are  obtained  for  the  dependent  variable  of  the  differential  equation.  Thus 
the  final  boundary  condition  cannot  be  fulfilled. 

The  linear  differential  equation  is  obtained  by  the  generalized  Newton- 
Raphson  method.  Therefore  these  unreasonable  values  of  the  parameters 
during  the  first  few  iterations  should  be  expected. 

In  order  to  avoid  these  unreasonable  values,  the  parameters  should  be 
changed  to  within  reasonable  ranges.  When  these  restrictions  are  used, 
the  convergence  problems  are  reduced. 

The  method  converged  to  the  exact  values  in  six  iterations.  Figure  (1) 
shows  the  nature  of  convergence  of  the  method. 


20 


PART  II 
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INTRODUCTION 

The  estimation  of  parameters  in  nonlinear  partial  differential  equations 
is  much  more  difficult  than  in  ordinary  differential  equations.  If  the 
model  equations  are  linear  an  analytical  solution  is  usually  obtainable 
from  which  the  parameters  may  be  estimated.  When  the  partial  differential 
equation  is  nonlinear  some  sort  of  extremely  time-consumming  trial  and  error 
integration  of  the  equations  becomes  necessary. 

Parameter  estimation  for  the  humidity  diffusion  in  concrete  has  been 
presented.  The  drying  of  concrete  is  one  of  the  basic  factors  in  creep, 
shrinkage  and  crack  formation.  In  the  construction  of  massive  structures, 
knowledge  of  the  rate  and  the  amount  of  heat  evolved  during  the  process  of 
hardening  are  desirable.  The  physical  process  is  assumed  to  be  represented 
by  a  nonlinear  parabolic  differential  equation. 

The  phenomenological  characteristics  of  the  drying  of  concrete  were 
studied  experimentally  by  Pihlajavaara  [1].  A  mathematical  equation  for  the 
moisture  dependence  of  the  diffusion  coefficient  was  proposed.  It  was  ob- 
served in  the  experiments  that  the  diffusion  coefficient  underwent  a  sig- 
nificant change  when  the  process  of  drying  proceeded  from  100  per  cent  to 
70  per  cent  relative  humidity  of  the  pores.  Definite  conclusions  have  not 
been  reached  so  far.  This  is  partly  due  to  the  lack  of  experimental  data 
for  the  drying  of  concrete. 

Bazant  and  Najjar  [2]  performed  numerical  analysis  with  the  diffusion 
coefficient  proposed  by  [1].  Trial  and  error  was  used.  Finite  differencing 
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of  the  space  and  the  time  derivatives  was  carried  through  based  on  the 
Crank  Ni col  son  method.  A  more  pertinent  form  of  the  diffusion  coefficient 
was  proposed.  It  contained  three  parameters. 

The  objective  of  the  present  work  is  to  estimate  the  parameters  in 
the  diffusion  equation  by  the  quasi! inearization  technique.  The  diffusion 
coefficient  suggested  by  [2 J  is  used.  Experimental  data  of  Abrams  and 
Orals  [3]  is  also  used. 

Integration  of  the  linearized  equation  is  carried  through  by  the  finite 
difference  method.  Both  the  space  and  the  time  derivatives  are  discretized. 
The  Crank-Nicolson  method  has  been  used  in  the  finite  differencing.  An 
algorithm  is  devised  based  on  the  least  squares  method.  Initial  guesses 
are  made  for  the  parameters  and  also  for  the  function  value  of  humidity 
as  a  function  of  space  and  time.  The  initial  values  chosen  for  the  purpose 
of  integration  comply  with  the  initial  boundary  condition  and  set  the  in- 
tegration constants  equal  to  the  unknown  parameters.  The  parameters  thus 
obtained  are  used  recursively  to  obtain  improved  values. 

The  error  in  the  experimental  data  of  [3]  is  considered  small  in 
magnitude.  The  mathematical  model  is  devised  to  fit  the  experimental  data. 
The  concrete  medium  is  assumed  isotropic  macroscopically.  It  has  been 
assumed  that  the  process  of  drying  takes  place  isothermally. 
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CHAPTER  I 

THE  BASIC  DIFFUSION  EQUATION 

Concrete  is  a  capillary  porous  colloidal  material.     It  has  a  wide  range 
of  pore  sizes.     The  complicated  physico-chemical  microstructure  allows  dif- 
ferent kinds  of  individual  flow  to  occur  simultaneously.     Vapor  diffusion, 
saturated  and  unsaturated  capillary  transfer  are  typical  examples  of  these 

f 1 ows . 

The  process  of  diffusion  is  a  result  of  random  particle  motions  within 
the  medium.  Matter  is  transported  from  one  part  of  the  system  to  another  as 
a  result  of  the  random  particle  motions.  The  flow  due  to  diffusion  is  not 
caused  by  external  forces  but  it  is  due  to  concentration  gradients. 

Inadequate  knowledge  of  the  actual  mechanisms  of  heat  and  mass  transfer 
during  the  drying  process  makes  the  formulation  of  a  detailed  kinetic- 
mathematical  model  difficult.  Generally  engineering  problems  that  have  a 
multiplicity  of  fundamental  laws  in  operation  are  nonlinear  mathematical 

models. 

Concentration  dependent  diffusion  in  an  isotropic  medium  is  governed 

by  the  Ficks  equation, 

!ld  .  J_  [C(Y  J  |Ij  0) 

at   ax  L  v  dy  3xJ 

or 

3Y 


3 


/  «  V(C(YJ  Wd)  (2) 
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Equation  (2)  is  expressed  in  vector  notation.  The  diffusion  coefficient 
undergoes  a  significant  change  when  the  drying  process  progresses  from  100 
per  cent  relative  humidity  of  the  pores  to  70  per  cent  relative  humidity  of 
the  pores. 

Pihlajavaara  [1]  proposed  the  following  mathematical  equation  for  the 
diffusion  coefficient: 


C[Yd]  =  C,  •  (1  +  a  Yj) 


(3) 


where , 

a  and  y  are  parameters, 

C-j  is  the  diffusion  coefficient  at  100  per  cent  relative  humidity  of 

the  pores, 
Y  .  is  the  relative  humidity  of  the  pores. 

Bazant  and  Najjar  [2]  proposed  the  following  equations: 


C£Yd]  =  C]  -  {1  +  (1  -  a)(l  -  Yd)Y} 


and 


C[Yd]  =  ^  •  a  + 


1 

-  a 

\ 

1  + 

r1-^ 

Y 

i 

,  l-h 

(4) 


(5) 


where, 

h  is  another  parameter. 

The  diffusion  coefficient  characterizes  the  diffusion  process  and 
represents  the  rate  of  drying.  Besides  its  dependence  on  the  moisture, 
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the  diffusion  coefficient  also  depends  on  the  properties  of  the  cement 
paste,  history  of  hydration,  moisture  concentration  and  the  ambient  temper- 
ature. 

Figure  (3)  illustrates  the  various  phases  of  drying  experienced  by 
concrete.  It  shows  a  plot  of  the  rate  of  drying  against  the  average 
moisture  content. 

The  medium  is  initially  at  saturation,  constant  temperature  and  con- 
stant external  condition.  AB  denotes  the  constant  rate  period.  The 
surface  of  the  medium  is  sufficiently  wet  during  this  phase.  It  enables 
to  simulate  a  free  water  surface.  The  surface  remains  wet  during  the 
period. 

BC  denotes  the  first- foil  owing  rate  period.  The  rate  of  evaporation 
from  the  saturated  surface  is  higher  than  the  rate  of  liquid  diffusion. 
The  surface  begins  to  dry  out.  The  effective  area  of  the  wetted  surface 
decreases  linearly  with  the  moisture  content  [17].  In  this  period  the 
rate  of  drying  is  approximately  a  linear  function  of  the  moisture  content. 
The  first  falling-rate  period  ends  when  dryness  is  reached  at  the  surface. 
The  moisture  content  will  be  in  equilibrium  with  the  moisture  in  the  en- 
vironment. 

CD  denotes  the  falling  rate  period.  It  is  characterized  by  sub-surface 
evaporation.  The  plane  of  evaporation  receedes  further  into  the  body. 
The  liquid  evaporates  from  a  relatively  dry  surface.  The  evaporation  rate 
depends  on  liquid  transfer  from  the  interior  of  the  mass.  The  liquid  trans- 
fer will  be  from  the  surface  into  the  environment. 
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The  rate  of  hydration  tends  to  retard  over  extended  time  intervals. 
This  is  because  of  the  dependence  of  the  diffusivity  on  the  pore  humidity. 
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CHAPTER  II 
LINEARIZATION  OF  THE  NONLINEAR  PARABOLIC  DIFFERENTIAL  EQUATION 

The  development  of  high  speed  digital  and  analog  computing  devices 
has  spurred  a  substantial  growth  of  the  mathematical  science  of  numerical 
analysis.  There  is  no  extant  theory  for  solving  nonlinear  partial  differ- 
ential equations.  Many  of  the  methods  used  are  suggested  from  procedures 
that  can  solve  linear  equations.  Extension  of  these  ideas  to  nonlinear 
partial  differential  equations  is  very  difficult  if  a  thorough  analysis  of 
stability,  convergence  and  error  is  carried  through. 

The  preliminary  step  involved  in  all  numerical  methods  is  to  intro- 
duce nondimensional  variables  in  place  of  dimensional  variables. 

f-fe^mf)  CD 

The  variables  Y,  t,  x  and  C[Y]  are  dimensionless 
where, 


1j     -    Y„n 

Y  =  5    en 


'^en 


, 


(td  -  y  ■  c, 
3 


x  ,,  t,,  Y  ,  are  dimensional  variables  of  space,  time  and  relative  humidity 
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of  pores  respectively.     Yen,  tQ,  Ld  and  C-j   are  dimensional   variables  of 
relative  humidity  of  the  environment,  curing  time,  half-length  of  slab 
and  the  diffusion  coefficient  at  100  per  cent  relative  humidity  of  pores 
respectively. 

The  dimensionless  diffusion  coefficient  C£Y]  is  given  by, 


where, 


C[YJ  =  a  +       (]   '  a) 

1   +  f3Y  (1-Y)Y 


1   -  Y. 


en 


is  a  parameter. 


(2) 


v  •      d        en 

1     "     Yon 

en 


Introducing  equation   (2)   in  the  basic  diffusion  equation  (1) 


£l- 

3t 


-  !_  (I    +     (1  -a)        }  H  ) 
3X   l\        1  +  byM-yW   9X  I 


(3) 


Thus, 


il «  jl  L  +    n-«) 

3t        3*  1 


il  +  /    +     (i  -  ex)   ■  1  £y_ 

1    +  BY(1-Y)YJ    9X        r       1    +  BY(l-Y)Yi    3X2 


(4) 


For  simplicity  in  the  derivation,  equation  (4)   is  written  as, 
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jY  ,.  3C  3Y       c  32Y 
3t   "    3X    3X  3)(2 


(5) 


Thus, 


2 
V  =  f£Y.  f£.  ■*-£»  a'     3'  y] 

3X 


(6) 


Equation  (5)  is  nonlinear  and  it  is  linearized  by  the  generalized 
Newton-Raphson  method  as  below: 
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c  +  +■  h 

where  the  subscripts  n+1  and  n  denote  the  (n+1)   iteration  and  the  n 


iteration  respectively. 


th  ,. 


Equation  (5)  may  be  written  for  the  n      iteration  as, 
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Introducing  equation  (8)  in  equation  (7), 
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Since, 


C  =  g(Y) 
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Introducing  equations  (13)  in  equation  (10) 
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On  cancellation  of  terms, 
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Equation  (15)  will  determine  the  particular  solution.  To  obtain  the 
homogeneous  solution,  the  equation  is 


Y '    - 

n+T  '  9Y 


n  '  9x 


.  11 
n  '  9x 


+  C  •  1?I 

n+1     "    9X2 


n+1 


h2c 

n 

f9Y^ 

2 

+   3C 
n        9a 

9Y9a 

VI 


[A 

W 

2 

+   9C 

9Y93 

n 

M) 

n       96 

92Y 


"    9X2  ^ 


3n+l 


^92C 
[9Y9Y 

n 

L3Xj 

2 

+   3C 

n       9y 

92Y 


"  .  9X2  ") 


'n+1 


fA 

n 

8X 

2 

9C   9^Y 


;\ 


n   9Y   ax2  nj 


n+1 


(16) 


2    2     2 
In  equations  (8)  to  (16)  the  equations  for  C,  f£  ,  ^-|  ,  |y^  ,  |yL 
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,  etc  are  as  follows: 
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Equations  (15)  and  (16)  are  rewritten  with  the  terms  all  rearranged  as 
follows, 
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Boundary  Conditions 

At  the  exposed  surface; 


Y(t)  =  Y   .(t)   . 
'    environment 


for  all  time  t. 
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At  the  center  of  the  slab. 
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for  all  time  t. 
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(i )  Initial  boundary  condition 

Y(0)  =  1,    at  t  =  0  for  all  x,  except  at  the  exposed 

surface.  (20) 

(ii )  Final  boundary  condition 


Y(tx)  =  b   at  t  =  t,.  . 
f  fina 


final 
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where. 


b  is  the  experimental  value  at  t  =  t£.     •,   . 

final 


Applying  the  generalized  Newton-Raphson  formula  to  equations  (18) 
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CHAPTER  III 
FINITE  DIFFERENCE  METHOD 

The  equations  (15a)  and  (16a)  of  the  preceeding  chapter  are  solved  by 
the  finite  difference  method.  Finite  differences  of  both  the  space  and  the 
time  derivatives  are  carried  through.  The  increments  in  Y  at  each  time  step 
can  be  computed  from  either  implict  or  explicit  methods.  The  explicit  method 
leads  to  a  numerically  unstable  solution  even  if  the  diffusivity  is  considered 
constant. 

Backward  or  central  differences  in  time  steps  reduce  the  instability 
of  the  problem.  For  both  the  schemes,  the  numerical  process  at  constant  dif- 
fusivity is  stable  [2].  They  also  allow  the  time-step  to  be  increased  or  de- 
creased as  desired.  The  dampening  of  error  in  subsequent  steps  is  stronger 
when  the  backward  differences  are  used.  The  central  differences  are  usually 
more  advantageous  because  of  their  accuracy. 

The  Crank-Nicolson  method  of  averaging  the  approximations  in  the  t 
and  (t+1)   rows  is  considered  more  accurate  and  has  therefore  been  used. 
The  explicit  and  implicit  methods  lead  to  discretization  errors  of 
0[At  +  (Ax)  J.  The  Crank-Nicolson  method  is  stable  for  all  values  of  the 

?  2.       2 

ratio  At/ (ax)  .  It  converges  with  a  discretization  error  of  0£At  +  (ax)  J  £18], 
Although  this  is  a  distinct  improvement  over  the  explicit  and  the  implicit 
methods,  the  computation  is  more  complicated  than  for  the  implicit  method. 

The  following  difference  equations  are  obtained  based  on  the  Crank- 
Nicolson  method, 
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Introducting  equations  (1)  to   (7)   in  equation   (15a)  of  the  preceeding 
chapter, 
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Yn(p+1 ,t+l )-2Yn(p,t+l )+Yn(p-l ,t+1 )+Yn(p-H  ,t)-2Yn(ptt)+Yn(p-1 ,t)^ 


2(AX)' 
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fac 

9a 


I  \   +   9C 

n  *    (an+ran)  +  Tt 


n  *   <W»n>  +  £ 


f  u3C 


;n+1(p,t+1)+Yn+1(P,t)       Yn(Ptt+1)+Yn(p,t)^ 
(  5  ?  ' 


(2) 


Writing  r  =  n 

(AX)2 


and  rearranging  the  terms  in  equation  (8)  above, 


Yn+1(p,t+1)   •  r  .  (i+  Cn  +  f  |n   •  I  •    {Yn(p,t+l)-Yn(p-l,t+l)+Yn(p,t) 


B*C 


1 


Yn(p-l,t)}  -  ±\    n  •  £  •    {Yn(p+l,t+l)   -  Yn(p,t+1)  +  Yn(p+l,t) 

d  T 


Yn(p,t)}' 


+  Yn+1(P+l,t+l)    .    r  .    (-1)   •    ]Cn  .  l  +  f 


n  *  4 


{Yn(p+l,t+l)  -  Yn(p.t+1)  +  Yn(p+l,t)   -  Yn(p,t)}] 


+  Y_+1(p-l.t+l)   •    r  •    (-1)    .   Cp   •  {-  Yn+1(p,t) 


+  \  ■    [Yn(p+1 ,t+l )-2Yn(p,t+l )+Yn(p-l ,t+l )+Yn(p+l ,t)-2Yn(p,t)+Yn(p-l  ,t)]. 


3C 

8a 


sC 


3C 


n   '    <Vran>  +lB     n   '    (^T^  +  F    "   <Yn+1-Yn)  + 


3C 

n   VTn+l~V    "    3Y 


n       2 


(Yn+1(p.t)-Y  (p.t+l)-Yn(p.t))]  +  \  •   [Yn(p+l,t+l)-Yn(p,t+l)+Yn(p+l,t)-YnCp,t)]' 
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A 


9Y9a 


re 


2 


(an+l"an}  +  fvt?  L  "    (Vl'Bn)  +  i 


r  w  3  c 

n  '   (VrV  +  ^2 


n      2 


(Yn+1(P.t)-Yn.(p,t+l)-Yn(p,t)) 


.  r      £C 
4  '   3Y 


IX+,(p+l,t+l)-Yn+,(p,t+l) 


n      L'n+1 


n+1 


+  Yn(p+l,t)  -  Yn(p,t)]   •   £Yn+1(p+l,t)   -  Yn+1(p,t)]  +  §  •   Cn 


£Yn+1(p+l,t)  -  2Yn+1(p,t)  +Yn+1(p-l,t)] 


n+1 


n+1 


(9) 


Equation  (9)  is  written  explicitly  below: 


Yn+1(p,t+l)   .  r 


-  +  a     + 
r        n 


1    -   a. 


l  +  e. 


Yn(p,t+l)+Yn(p,t)1Yn 


V1 

Tn  r       Yn(p,t+1)+Yn(P,t)- 

|.  (1"a"h"&"    i1  "     ^- — -Hr--{Yn(P,t+i)-Yn(p-l,t+i)+Yn(P,t) 


Y. 


1  + 


1 


Yn.(p.t+1)+Yn(p.t)^2 


1       (1-a  h  & 
-  Y  (p-l.t)}"i 


n 


8 


y  -2 
Yn(p,t+1)+Yn(ptt)^  n 


Y. 


1    + 


Yn(p,t+1)  +  Yn(p,th   n 
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Y„  f  Y„(P.t+l)     +Yn(p,t) 


I 


(VD-  (vr) 


l:+e"  •    1 


Yn  3 
Yn(p,t+1)  +  Yn(p,t)^  n! 


Yn(p+l,t+l)   -  Yn(p,t+1)  +  Yn(p+l,t)   -  Yn(p,t)}2j 


+  Yn+1(p+l,t+1)    •   r   •    (-1)   •     \an  + 


1 


1   +  3. 


1 

YJ       2 
Y„(p,t+1)  +  Yn(p,thn 


+  i.(hA 


1 


Yn(P»t+D     +    Y     (P,t)v    n 


\7 

,         Yn    r         Y   (p,t+l)  +  Y   (p,t)1 


•   {Yn(p+l,t+l)  -  Yn(p,t+D  +  Yn(p+l,t)   -  Yn(p,t) 


+Yn+1(P+l,t+l)    •   r   •    (-1)    •   L    + 


1-c 


1    + 


1    - 


Yn(p,t+1)  +-Yn(p,th   n 


Yn+1(p,t)  +  £  •   £Yn+1(p+l,t+l)  -  2Yn(p,t+l)  +  Yn(p-l.t+l)  +  Yn(p+l,t) 


2Y 


n(p.t)  +  Yn(p-l,t)] 
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l-> 


T„    (        Yn(p,t+1)  +  Yn(p.th   " 


1   + 


(a    ,,    -   a„) 
v    n+1  IV 


n-«n)Vnn 


-1 


Y„(p.t+1)  +  Yn(p.th   " 


W( 


Yn(p,t+1)  +  Yn(p,t). 


n  t 


<Vl   "  3n> 


o-A" 


Yn(p,t+1)  +  Yn(p,t)1   n 

•   in 


1   - 


Yn(p.t+1)  +  Yn(p,t) 


r        y. 

41  +  e_r 


Yn(p,t+1)  +  Yn(p,t)^Tn 


<Vl   "  V  + 


1 
2 


yn  r       Yn(P,w)  +  Yn(P.t)1 


V1 


n'  'rrn 


1   + 


1   - 


Yn(p,t+1)  +  Yn(p,t)< 


rrtW^.-V**1)  -V**?) 


+  £■   £Yn.(p+l,t+l)  -  Yn(p,t+1)  +  Yn(p+l,t)   -  Yn(p,t)r 
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-Y    B 

rrrn 


f        Yn(p.t+1)  +  Yn(p,t)1   n 


Y„-l 


{Vi  -  %] 


1 :+  a. 


1  - 


Yn,(p,t+1).+  Yn(p.t)^   n 


o-n^V 


Y.-l 


1    - 


V1 

Yn(p.t+1)  +  YJp.th 


1  ■  + 


Yn(p,t+1)  +  Yn(p,t) 


Y     3 

1  n 


i  -  e. 


Yn(p,t+1)  +  Yn(p,t)/n 


(e«*i  -  e«)  + 


i  +  e. 


Y     3       VPn+l       *V 
Yn(p,t+1)  +  Yn(pst),   n> 


H»C   ' 


Yn(P.t+U  +  Yn(p.t) 


1   +  B. 


1    - 


Yn(p,t+1)  +  .Yn(p,th   " 
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Y  (p,t+l)  +  Y(p,thYn 


.  {I  -  Y     •  In 


Yn(p,t+1)   +  Y   (p,t)^        t 
1   -J -2 2 W     l+yn   •   in 


Yn(p.t+1)  +  Yn(p,th^ 


j    '    ("Vl   "  Yn} 


v       ir'n      pn 


f        Y(p,t+1)  +  Y_(p.th   n 


1 


n 


1   + 


T     3 
Yn(p,t+1)  +  Yn(p,t)^   n. 

2  J     J 


h  -  n 


Y„(p.t+1)   +  Yn(p,th    n 


2 


(YB+T) 


(Y„-D 


1    + 


n     '    I 


Yn     3 
Yn(p,t+1)  +  Yn(p.th 
2 


I'  {Vl^'^-VP^^  *  Yn(p't}} 
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+  i 


V1 

Y  ,    Y  (p.t+1)  +  Y  (p.t)-i 

(i-Oy„Oi  -J ?-* 


n'  'n  n 


Tn  r    Yn(p,t+1)  +Yn(p,t)/n 


T 


[Yn(p+l,t+1)  -  Yn(p,t+1)  +  Yn(p+l,t)  -  Yn(p,t)] 


CW1*1^  "  Vl(P't)] 


n-  an  + 


1  -  a. 


1  + 


Yn(p,t+1)  +  Yn(p,t)/n 


J 


CYn+1(P+l,t)  -  2Yn+1(p,t)  +Yn+1(p-l.t)] 


n+1 


n+1 


(10) 


If  the  space  domain  is  divided  into  M-l  equidistant  spaces  by  M 
nodes,  there  will  be  M  equations  corresponding  to  equation  (10).  The  first 
and  the  last  of  these  equations  satisfy  the  boundary  conditions, 


W^'^enviS^ 


ronment 


(Ha) 


and 


3Yn+1(M,t) 


3X 


=  0  - 


Thus 
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2ax 


and 


W*1-^  =  Vi^1^ 


(lib) 


The  system  of  M  simultaneous  algebraic  equations  is  written  in  matrix 
form  as, 


[A]Mvm  '    mMvl   =  [K] 


JMxM 


Mxl        L'NJMxl 


(12) 


where, 


[A]..  M  is  a  tri diagonal   band  matrix, 


fYl.  ,    is  the  column  matrix  with  relative  humdiity  predictions 
at  the  subsequent  time  step, 

[K]..  ,   is  the  column  matrix  with  elements  of  known  values. 

Typical   values  that  constitute  the  p      row  of  the  matrices  [A]  and 
[K]  are  written  below 


A[p,p-2]  =  0, 


p  =  1  ,M 


(13a) 


A[p,p-1]  =  r  (-  1)    •  L     + 


1    -   a. 


1    + 


Yn(p,t+1)  +  Yn(p,t)/n 


P=l  ,M 


(13b) 
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A[p,p]  ■  r  •    ]f  +  ^n  + 


1    -  a_ 


1   + 


UP.t+1)  +  Yfp.t) 


1 
4 


(1-a    )y   B 


r         r         Yn(p,t+1)  +  Y^p.thV1 


1   +  6. 


Yn(p,t+1)  +  Yn(p,t) 


Yn.2 


{Yn(p,t+1)   -  Yn(p-l,t+l)  +  Yn(p,t)   -  Yn(p-l,t)} 


I 

8 


H)W"   •    f1    -Yn(P.^)^Yn(p,t)]V2 


1    +  6. 


1    - 


Yn(p.t+1)  +  Yn(p.t) 


Y„.3 


Yn(p.t+1)  +  Yn(p,t)/n 


(VD-  (y„-D 


i  +  e 


n 


1  - 


Y„(p,t+1)  +  Yn(p,t)/n^3 


|Yn(p+l,t+l)   -  Yn(p,t+1)  +  Yn(p+1,t)   -  Yn(p,t)}2 


,     P  =  1,M 


(13c) 
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A(p,p+1)  =  r  •    (-1) 


1    -   a 


an  + 


1  + 


Up.t+1)   +  Yn(p,t) 


1 
2 


x       n     n  n 


1   - 


Y„(p,t+1 )  +  Yn(p,t)1   n 


Y„-l 


1   + 


1   - 


Yn(p,t+1)  +  Yfp.t)/^2 


1 

4 


{Yn(p+l,t+l) 


Yn(p,t+1)  +  Yn(p+l,t) 


Yn(p,t)} 


p  =  1,  M  (13d) 


A(p,p+2)  =  0, 


p  =  1,  M 


(13«) 


K[p]  -  Yn+1(p,t)   +§  •    [Yn(p+l,t+l)   -   2Yn(p,t+l)   +  Yn(p-l,t+1)   +  Yn(p+l,t) 


2Yn(p.t)  +  Y  (p-l,t) 


1   - 


1 


yn       r        Yn(p,t+1)  +  Y  (p,t) 


r}  •  <vi  -  «n> 


"-nVn 


Yn(p,t+1)  +  Yn(p,t)/n 


1  + 


Yn(p,t+1)  +  Yn(p,t)Jn^2 


<pn+l    "  *n> 
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U-nK 


Yn(p,t+1)  +  Yn(p,t) 


in 


n 


Yn(p.t+1)  +  Yn(p.t)- 


1   + 


Yn(p,t+i)  +  Yn(p,th"V 


<Vl    "   Yn)    + 


(l-o    )y  B        • 
v       n/rrrn 


Y„(p,t+1)  +  Yn(p,thV1 


Y. 


enn  •  h 


Y„(p,t+1)  +  Y   (p,t)/n>2 


n 


1 
2 


(Yn+1(p,t)   -  Yn(p,t+1) 


Yn(p,t)} 


[Yn(p+l,t+l)   -  Yn(p-.t+l) 


+  Yn(p+l,t)   -  Yn(p,t)]' 


Y. 


nun 

'npn         | 


,        Y  (p,t+l)  +  Yn(p,t)Jn"1 


1  +  6. 


1   - 


Y„(p,t+l)  +  Y(p,thV   '    (%+1   "^ 


o-rf/1 


Yn(p,t+1)  +  Yn(p,t)/n 


y  -1 


Y. 


1    + 


Yn(p.t+1)  +  Yn(p,thTm3 


J     J 
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1   - 


n 


1    - 


Yn(p,t+1)  +  Yn(p.t) 


n 

Li 


Y         {        Yn(p,t+1)  +  Y   (p,t)/n 
1  +  *n     '    P  '  ~ 2 — 


yr  3 


<Vl   "  en>  + 


v         f        Vn(P,t+1)+Yn(p,tyn 


Y.-l 


{'  ♦  -*" 


1    - 


Yn(p,t+1)  +  Yn(p,t)^ 


r,  3 


1  - 


Y„(p,t+1)  +  Yn(p,t) 


1   -  Yn   •   in 


V     ] 


Yn(p,t+1)  +  Yn(p,t)V 


+  \l   +  y_  In 


n       I 


r       Yp(P>t+1)  +  Yn(p,t) 


(Vl  "  Yn) 


(1-a   )    •   Y  6 
v       n'       'n  n 


1 


Yn(p»t+1)  +  Yn(p,t)Jn"2 


2 


1   + 


y  Y    (p.tfrl)    +    Y    (P,th  V 

3  "    •      1 
n 


Y„(p,t+1)  +  Y   (p,t) 


(Yn+D    -<Yn-l) 


1    + 


Y  (p.t+1)  +  Yn(p,t)^n>: 
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1  •    {Yn+1(p,t)   -  Yn(p,t+1)  -  Yn(p,t)} 


r      "'%K  '  *„n 


1    - 


Yn(p.t+1)  +  Yjp.thV1 


Yn        f         Y   (p,t+l)  +  Y   (p,th"S2 
1  +  enn  •     1 


CYn(p+l.t+l)  -  Yn(p,t+1)  +  Yn(p+l,t)  -  Yn(p.t)]   •    [Yn+1(p+l,t)  -  Yn+1(p,t)] 


n+1 


n+1 


♦* 


1    -   a. 


a     + 

n 


1  + 


n         ( 


1   - 


Yn(p,t+1)  +  Yn(p,t) 


[Yn+1(p+l,t)  -  2Yn+1(p,t)  +  Yn+1(p-l,t)] 


(14) 


Matrix  [Y]M  -,   is  composed  of, 
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fW1^1)  1 


Yn+1(2.t+l) 


Yn+1(p,t+l) 


Yn+1(M-l,t+l) 

lWM't+1)  J 


(15) 


Mxl 


The  initial  values  for  obtaining  the  particular  solution  are, 


W1-01 


YP,n+l(2'0) 


Yp,n+l(3''°) 


W^'^ 


Vn+l^1'*) 


ap,n+l(0) 


p,n+l 


rp,n+l 


(0) 
(0) 


(y      .(o)      A 
environment 


1 

0 
0 
0 


(16) 
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Equation  (16a)  of  the  preceeding  chapter  is  used  to  determine  the 
homogeneous  solution.     It  is  rewritten  below 


wi       _  aC 

n+1  "  dy 


3Y    I      <    aY 

n  '   ax  'n       ax 


+  C. 


a2Y 


n+1        n      ax2 


n+1 


aY 


SX 


ra2C 


aYaa 


n   "   an+l       3Y36 


2 

n       Bn+1       3Yay 


a2C 
n  '  Vl  +  ^2 


n  *    'n+1 


ax2 


aC 

3a 


+   9C 

n        n+1       ae 


n       pn+l        ay 


Y 


aC 


n       rn+l       aY 


n  '     n+1 


(17) 


Using  the  finite  differences   in  equations   (1)  to  (7)  in  equation  (17) 
and  proceeding  in  a  similar  manner  as   in  equation  (10),   the  final   form 
obtained  is, 


Yn+1(p,t+l)   .   r  . 


1    -  a. 


-  +  <a     + 
r       in 


1  + 


1   - 


Yn(p,t+1)  +  Yn(p,t)/n 


+  -f 


(1-a    )y   6   "    • 
1       v       n;rn  n 


1 


Yn(p,t+1)   +  Yn(p,thV] 


Yn        t         Y   (p,t+l)  +  Y   (p,t)/n^2 
1  +  C   •    V   -  ~ 2 


[Yn(p,t+1)   -  Yn(p-l,t+l)   +  Yn(p,t)   -  Yn(p-l,t)] 
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1        H-an)YnBnn 


Y  (p,t+l)  +  Yn(p,th 


Y.-2 


1   +  P. 


1    - 


Yn(p,t+1)   +  Y^p.t)^ 


T>w)*Vp.tV-.(     ,.(    „ 


r.      y, 


i1  +  v-  ] 


Yn(p,t+1)  +  Yn(p,t)/ni 


{Yn(p+l,t+l)  -  Yn(p.t+1)  +  Yn(p+l,t)  -  Yn(p,t)}2] 


+  Y  +1(p+l,t+l)    •   r  .'  (-1)    .       an  + 


1    -    a. 


Y. 


'   ♦  »„"   •    I'   -  ^ 2 


1 

Yn(p.t+1)   +  Yn(p,t)V 


(l-an)ynBn     •      1    -  - 


Y„(p.t+1)  +  Yn(p,t)/n"' 


* 


1    +  B.."    •      1    -      " 


Y„(p,t+1)  +  Y   (p,t)/n^2 


J  •    {Yn(p+l,t+l)   -  Yn(p,t+1)   +  Yn(p+l,t)   -  Yn(p,t)}j 


+  Yn+1(p-l,t+l)   •   r  .    (-  £) 


1    -    a_ 


an  + 


1   + 


Yn(p,t+1)   +  Yn(p,t)/n 
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Yn+1(P,t)+f 


Yn(p+1,t+l)   -  2Yn(p,t+l)  +  Yn(p-l,t+l)  +  Yn(p+l,t) 


2Yn(p,t)  +  Yn(p-l,t)] 


1 


1   + 


YJp.t+l)  +  Yn(p,t) 


an+l 


Y  -1                 Y   (p,t+l)  +  Yn(p,t)  n 
(1-ajY.e"        •      1-J 5 ~ 


n"n  n 


1  + 


Y„(p,t+1)  +  Yn(p,t)Jn12 


>n+l 


Yn       r        Y  (p,t+l)  +Y  (p.t)/" 
(l-an)enn  •    II  -ir-  •   An 


n 


1 


Yn(p,t+1)   +  Yn(p,t)- 


1   + 


\  2 
Yn(p.t+1)  +Yn(p,th   "> 
_ 


Vl  + 


(1-a   )y   6  n 


1   - 


Yn(P.t+l)  +  Yn(p,t)/n"1 


1    +  6. 


Yn(p,t+1)  +  Yn(p,thYrh2 


{'   Yn+1(p.t) 


+  £  •    [Yn(p+l,t+l)   -  Yn(p,t+1)   +  Yn(p+l,t)   -  Yn(p,t)]' 
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n  n 


7  -1 
Yn(p.t+1)  +  Yn(p,t)i 


n 


1  + 


2  J     J 


r  2  '  an+l 


0-a    )y26  n 
v       ir'n 


Y.-l 


Yn(p,t+1)  +  Yn(p,thV 


1  + 


Y. 


Yn(p,t+1)  +  Yn(p,t)^ 


yr  3 


1    -   B. 


1 


Yn(p,t+1)  +  Yn(p.t) 


Y          ,         Yn(p.t+1)   +  Y  n(p,t)/ni3 
!  +  3n"   •    (l   -  J H 


"-A" 


n       f,       Yn(p't+1)  +  Yn(p>t)Vn 

2 


1  + 


Y--1 


Yn(p,t+1)  +  Yn(p,t)i 


V 


Yn(p,t+1)   +  Y(p,t)- 


[1   -  Tn  •   in{gn 


f         Y„(p,t+1)   +  Y(p,t) 


1    - 


2 
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^n  *n  1*n  ' 


Yn(p.t+1)  +  Yn(p,t)1 


Vl  + 


v      n/Tn  n 


1  - 


Yn(p,t+1)  +  Yn(p,t)/n" 


Y  Y   (p,t+l)  +  Y   (p.th 

1   +  Bnn   .     !    .  J] _i! 


Y-  3 


1    - 


Y„(P.t+1)   +  Yn(p,t) 


(VD-(Yn-D    •!•  Yn+1(P,t) 


(1-a   )   •   Y  e 
r       v       n;       rn  n 


1 


Yn(p,t+1)  +  Yn(p.t) 


Y„-1 


r  Y„        r         Yn(ptt^)+Yn(p.t)Jn,2 

I1  +  *n     *    i1 ^ J 


[Yn(p+l,t+l)   -  Yn(pft+1)  +  Yn(p+l,t)   -  Yn(p,t)]    •    [Yn+1(p+l,t)   -  Yn+1(p,t)] 


n+1 


^ 


1    -    a_ 


*n  + 


1  + 


Yn(p.t+1)  +  Yn(p,t)/n 


[Yn+1(p+l,t)   ■  2Yn+1(p,t)  +Yn.1(p-l.t)] 


'n+1 


n+1 


n+1 


(18) 
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Equation   (18)  also  consists  of  a  system  of  M  simultaneous  equations. 
In  matrix  form  the  simultaneous  equations  are  written  as, 


^MxM  '    ^Mxl   ~   ^'^Mxl 


(19) 


J.L. 

The  elements  in  the  p      row  of  [A]  and  [Y]  are  identical   to  those 
in  equation  [12].     However  the  elements  in  the  p       row  of  [K']  are  different 


K'[p]  =  Yn+1(p,t)  +  £  •    [Yn(p+l,t+l)   -  2Yn(p,t+l)  +  Yn(p-l,t+l)  +  Yn(p+l,t) 


■2Yn(p,t)  +  Yn(p-l,t)]   • 


1   + 


Yn(p.t+1)  +  Yn(p.t)/nJ     n+1 


(1-a  h  e  n 
v       n;tn  n 


1   - 


Yn.(p,t+1)  +  Yn(p,t)/n 


1  + 


Yn(p,t+1)  +  Yn(p,t)^ 


Y,   2 


>n+l 


V'%K 


Yn(p,t+1)  +  Y  (p,thYn 

■J = 2 .   in 


/        Yn(p't+1)  +  Yn(p>t)^ 


Yn(p,t+1)  +  Yn(p,t)1   n! 


rn+l 


1   + 


; 
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<W   •  W 


1      - 


Y„(p,t+1)  +  Yn(p,t)/n-1 


1   +  6. 


Yn(p.t+1)  +  Yn(p,t+l)Jm2 


i'Yn+1(p.t) 


+  J-  •    CYn(p+l,t+l)  -  Yn(p,t+1)  +  Yn(p+l,t)   -  Yn(p,t)]' 


Yn(p,t+1)  +  Yfp.tjJrf1 


r   y 


n  n 


1  + 


Yn(p,t+1)  +  Yn(p,t)/n 


7  •  Vl  + 


(l-an)7n8 


2  Vl      f        Yn(p,U1)  H^tl/n-1 


(1  + 


Yn(p,t+1)  +  Y^p.thV3 


r        Yn(p,t+1)  +  Yn(p,thYn 


Y, 


1    +B*     1    - 

n 


Yn(P.t+l)   +  Yn(p,thV  n+1 
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»-»>»  " 


1   " 


Y„(p,t+1)  +  Yn(p5t)Jn_1 


1  +  3 


n 


Yn(p,t+1)  +  Y  (p,t)Jn^3 


Yn  r         Yn(P>t+1)  +  Yn(p,t)/n 
Bn  '  '  '      2 


l-Yn«HPn 


1  - 


Yn(p,t+1)  +  Yn(p,t)^^ 


+  Yn  •  *n  jgn  .|1 


Y  (p,t+l)  +  Yn(p,t)v 


rn+l 


n    n   n 


Yn(p,t+1)  +  Yfp.thV2 


1  +  6. 


n  n 


Y  >,t+l)  +  Y  „(p.thYih3 


f    Yn(p,t+1)  +  YnCp,t) 


1  +  B. 


r   y  (p,t+i)  +  y  (p.t)/^3 


'  •  (Y  +1)  -  (Yn-l)l    ,  1 

D . D L  .  1  .  y   (Dt) 


2   n+1 


(20) 


1 


There  are  three  sets  of  homogeneous  solutions.  The  initial  values  for 
obtaining  these  are  given  below, 
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Wl(1'0)  Wl(1'0)  Wl(1'0) 


YhUn+l<2'0)  Yh2W2>0>  Yh3,n+l(2'0) 


Yhl,n+l(3'°)  Yh2,n+l(3'°)  Wl<3'0) 


Yhl,n+l(P>°>  Yn2WP'0)  Yh3,n+l(p*0) 


Wl(,M*0)      Wl(W'0)      Wl^1'0* 


hl,n+l 


ahl,n+l 


5hl,n+l 


(M.0)  Yh2     +1(M-0) 


(0) 


(0) 


ah2,n+l(0) 


Wl("'0) 


tth3,n+r(0) 


Wl(0)  eh3,n+l(0) 


*hl.n+l<0)  ^h2,n+l(°)  ^h3,n+l(0) 


0       0       0 


0       0       0 


0       0       0 


0       0       0 


c 

0 

0 

0 

0 

0 

1 

0 

0 

0 

1 

0 

0 

0 

1 

(21) 


The  general   solution  of  the  system  of  equations   is  obtained  by  the 
principle  of  superposition 


W^  =  Vn+l"'**  +  kl}   ak,n+l    •   Wl"'**'     j  "  ] 'M 


(22a) 


an+1(j,t)  =  o__n+1(t)  +     I.   a 


p.n+P1"       |c41  ak,n+l       uhk,n+l 


(t) 


(22b) 
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WJ'tJ  "  W4)  +  ,£  ak,n+l    *   •Hik.nfl^)  (22C) 


WJ,t)  =  Yp,n+l(t)  +  kl}   ak,n+l    •   Yhk,n+1 


(t)  (22d) 


The  subscripts  p,  hi,  h2  and  h3  denote  the  particular  and  three  sets 
of  homogeneous  solutions.     The  initial   values  for  obtaining  these  four  sets 
have  been  selected  to  satisfy  the  initial  boundary  condition  and  also  to  set 
the  integration  constants  equal   to  the  parameters. 

It  may  be  seen  from  equations  (22b),    (22c)  and  (22d)  that, 

an+1(0)  =  a1>n+1  (23a) 

W»  =a2,n+l 
and 

Vi'°l  =  Vi  (23c) 

Since  the  parameters  are  assumed  constant  functions,  therefore, 


(23b) 


an+l(t)  =  al,n+l 

en+1(t)  =  a2>n+1 


(24a) 
(24b) 


vi^-Vi  (24c) 
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CHAPTER  IV 
NUMERICAL  ANALYSIS 

The  linear  differential   equation  is  solved  in  two  steps,  assuming  that 
a  solution  exists  for  the  problem.     First  one  set  of  particular  and  three 
sets  of  homogeneous  solutions  are  obtained  numerically.     Initial   values  for 
obtaining  each  of  these  sets  have  been  stated  in  the  preceeding  chapter. 
The  tridiagonal  band  matrix  is  inverted  at  every  time  increment  by  using  the 
Thomas'  algorithm.     This  algorithm  expedites  the  numerical  work  involved  at 
each  time-step. 

The  next  step  involves  the  determination  of  the  three  integration  con- 
stants.    The  least  squares  method  is  used.     The  error  between  the  experi- 
mental  data  and  the  model   prediction  is  minimized. 


m,     M 
1     M  o 


ml     M 
Qn+1  =  ^  JJ1^n*l«-V+al,n*l   -Wl«-V 


+  a2,n+1    •   Vz.n+l^V  +  a3,n+l    '   Yh3,n+l(J'V   "  W^S^ 

(2) 

Since  all   the  particular  and  homogeneous  solutions  are  known, 
equation  (2)   is  written  as, 
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ml     M 
Qn+1  °  sl}    J}  ^.n+l^V  +  al,n+l   '   ^.n+l^'V 

+  a2,n+l    '   ci3,n+l(j'ts)  +  a3,n+l    "  ^.n+l^V^  (3) 

Minimizing  equation  (3)  wrt  a-.     +1 ,  a2    +1   and  a3  n+1   the  following 
three  equations  are  obtained, 

ml     M 
Ji    .^  q2.iH-l(J,ts)   '   Cql,n+l(j'ts)  +al,n+l    '  q2,n+l(j>ts} 

+  ^l^'V^n^'V  +  a3,n+l    ■  Vn+l"'*^  =  °        (4a) 


m,     M 
sl}  JS1  q3.n+l(J'ts)   '  hl.n+lU'V  +  al  ,n+l   *  ^.n+l^'V 


+  a2,n+l   '   q3,n+l^V  +  a3,n+l    '   q4,n+l  (j'VJ  =  °  (4b) 


ml     H 


l}    -l}  'U.n+lU'V  '  *h.n*l{J'ts)  +  al,n+l   '  q2,n+lCj»V 


s=i  J 


+  a2,n+l    *  "s.n+l^V  +  a3,n+l    '  q4,n+l^VJ  =  °  (4c) 
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The  three  simultaneous  algebraic  equations  (4)  are  solved  to  obtain 
integration  constants.  They  are  equal  to  the  unknown  parameters.  The 
parameters  are  used  recursively  to  obtain  improved  results. 

The  flow  chart  figure  (4)  that  was  used  in  the  numerical  solution 
appears  in  the  Appendix. 
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CHAPTER  V 
CONCLUSION 

In  the  solution  of  the  problem  convergence  was  not  obtained.  The 
following  factors  may  be  attributed  to  this  behavior. 

Actual  experience  has  shown  that  the  generalized  Newton-Raphson  method 
is  very  sensitive  to  the  error  of  the  unknown  or  guessed  initial  condition, 
and  is  unstable.  The  guessed  values  of  the  unknown  parameters  as  well  as 
the  functional  equation  of  humidity  for  the  first  iteration  have  to  be  very 
close  to  the  correct  values  for  the  problem  to  converge. 

The  lack  of  convergence  encountered  in  solving  the  problem  may  also 
have  been  due  to  insufficient  spacial  nodes.  Numerical  analysis  with  upto 
thirty-four  nodal  points  in  the  spacial  coordinate  and  one  hundred  in  the 
time  coordinate  was  carried  through. 

The  Newton-Raphson  method  can  converge  quadratically  only  if  the 
method  should  at  all  converge  [11]. 
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l.Q 


s  0.5 


o 

O 

o 


0.4 


0.3 


0.2 


0.1 


0 


THIRD  ITERATION 


_i I L 


J 1 I I L 


0    0.1   0.2   0.3    0.4   0.5   0.6   0.7   0.8    0.9   1.0 

AXIAL  LENGTH,  t 

FIGURE  1.  CONCENTRATION  PROFILES  WITH  AXIAL  DIFFUSION, 

SHOWING  THE  NATURE  OF  CONVERGENCE  OF  THE  METHOD. 
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START 

t 

READdx,dt,N,E 

t 

READ  Yn+I(t,i),t=0,T;i»I,M 

PARTICULAR  SOLUTION 
Vlf0.1)  "  I!   1»2.«J  Vl=0;   Sn+I=0i  Vl"0;  W°'1J 


SOLVE  THE  EQUATIONS  BY    jjSIUG  ^  ?otffitK  0RDB 
RTJNGE-KUTT4  GILL  METHOD 


N=I 


N-2 


Yp(t.D=Yn+I(t,i) 
t-0,T  1«I,M 


t=0,T     i=I,M 


HOMOGENEOUS  SOLN  I 
Vn+I(0,i)-0;1-I,M 

N«2 

3 


HOMOGENEOUS  SOLN  II 
Yn+I(0,i)-O;i=I,M 

■*r0'*n+rI»TB+i"( 

N=3 


t*0,T  i-I.M 


HOMOGENEOUS  SOLN  III 
N=4 


t-0,T;  1-I.H 


SOLVE  THE  LEAST  SQUARES  EQUATIONS  AND  OBTAIN 
Vl'Vl  AND  Vl 


I 


FORM  THE  GENERAL  SOLUTIONS 

Vl  <*-1  ^p.n+I  <».*  >+VlYhl ,n+I  ^  >+*n+IYh2,n+I  tW  >*WW,  U.  1 ) 
t«0,T;  1-1.M 


>£ 


lVi<*»1>*YniiTA(W)|;  « 


T 


<£ 


t-0,T  1-1, M 


ISTOP 


FIGURE  2.   FLOW  CHART  FOR  SOLVING  THE  NONLINEAR 
ORDINARY  DIFFERENTIAL  EQUATION. 
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A 

b\ 

i— 

■zi 

LU 

1— 

^y 

o    - 

o 

UJ 

CI 

ZD 

1— 

oo 

H- 1 

o 

UJ 

cV 

CD 

<C 

c2 

UJ 

- D 

> 

< 

RATE  OF  DRYING  — »- 

-« —  MOISTURE 

FIGURE  3:  PLOT  OF  THE  RATE  OF  DRYING  AGAINST  THE  AVERAGE 

MOISTURE  CONTENT. 


START 
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READdx.dt.N.e 


I 


READ 


Yn+I(t,i),t=0,T;i«I,M  J 


PARTICULAR  SOLUTION 
Yn+I(0,i)  -  1;   1-2.H;  an+I»0;  Vl»0;  rn+I»0;  Y„+I<0.1) 


SOLVE  THE  EQUATIONS  BY   INVERTING  THE  TRIUIAGONAL 
BAND  MATfelX 


N-I 


N*2 


t*0,T  1«I,H 


t-0,T     i=I,M 


t 


HOMOGENEOUS  SOLN  I 
Yn+I(0,i)=0;i»I,M 

Vl^'Vl^'Vl"0 
N"2 


HOMOGENEOUS  SOLN   II 
Yn+I(0,1)»0;1-I,M 


Yh2(t,i)-Yn+I(t,i) 
t*0,T  1-1, M 


ln+I 


5n+I 
N=3 


HOMOGENEOUS  SOLN  III 
Yn+I(0,i)-0;i«I,M 

'an+Ix°'sn+r°'VlBl 
N«4 


-I  r— 


yh3(t^CVl(t,1) 
t»0,T;  1»I,M 
I 


SOLVE  THE  LEAST  SQUARES  EQUATIONS  AND  OBTAIN 

Vi'Vl  AND  Vi 


I 


FORM  THE  GENERAL  SOLUTIONS 
t«0,T;  1-l.M 


iVl(t'1')-YDATA(t-i)l;   £KJ  Yn(t.1)-»n*i(t.iJ 

t«0,T  1»I,M 


lJ 


STOP 


FIGURE  H.      FLOW  CHART  FOR  SOLVING  THE  NONLINEAR 
PARABOLIC  DIFFERENTIAL  EQUATION. 
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LEVEL  21.7    «    JAN   73    I  CS/360      FORTRAh  H  BATE      78. 1 15/22.21.37 


THIS    PROGRAM    ESTIMATES    PARAMETERS     IN    A    NOIL  INEAR    ORDINARY 
DIFFERENTIAL    EQUATION.     THE    CENEK4LISE0    KEWTCN-RAPMSCN    XtTHOO 
IS    USED   TO    LINEARISE    THE    fc'.'UATJUN.    USt    L>    THE    RUNGE-KUTIA    CILL 
METHOD    IS   MACE    FOR    NUMERICAL    INTEGRATION. 

DIMENSION  XO(101),YOtlOl),ROtl),PO(ll,XntCI)  ,TlllCUfRI(l»  fPltl), 
•XPl 101) ,YPC 101 ) ,XH 1(10 H  .XH2U01  ),XH3(101>  ,  Yhl  ( 1 01 1 ,  YH2  U3  11 .  YH3  U 
•  01  >  tttSt  1G)  ,  Wl  <■)  ,'r  I  28)  ,Cn  Ctl  ,A[  3 J 

COKKCN/  F  Cftfcfi  R/  X  C ,  Y  0 ,  RO ,  PO 

COMMON/lCUNTk/INT.ICHECK 

1  FORMAT  IF5.2,2X,F5.l,2X,I4,2X,I2.2X,  I  3,  2X  ,F  10.  fc,2X,F8.51 

2  FORHATIF8.5) 
J   FORMAT  I/// J 

A    FORMAT  m.'RO1,    5X, F3 . 1 , /, 9X, 'PO ' ,     5X, F3 .1 , // , T12, 'N' , T2S , 'R« ,T32, 

••P' ,TA3, 'GENERAL    SOLUTIONS    AT    CORRESPCNCINC    DATA    POINTS', /I 
5    FORMAT  I IH1) 
13    FORMAT (2F10. 6) 

63    FORMAT! '0'  ,£X, I 3. 5X.2F9. 5, SX , 10F5.21 
WRITEI6.5) 

READ (5,1!    DT.TF.M, N.NDATA.XE.EPS 
K1«M*1 

DO    10    1-l.NOATA 
READI5.2)    BSUJ 
10    CONTINUE 

97   CONTINUE 

READ(5,13)    RINT.PINT 

KA«1 
9 A    CONTINUE 

REACIS.13)    A2.A3 

VRITE(6,<i)    RINT.PINT 

DO    12    l-l.Xl 

X31 11*0.83129 

YOt  11—0.5 
12    CCNT INUE 

R0tl)«RINT 

P0(1)*PINT 

NN=1 
90   CONTINUE 

xim«xc 

YKU-O. 
Rltll-O. 
Pit l)«0. 
IDUPKY«1 
ICHECK'l 
20   CONTINUE 
U(l)«Xl(ll 
WI21-YK  11 
W131-R1UJ 
Vl4)«Pltl) 

c 

C  INTEGRATE    THE    DIFFERENTIAL    EQUATION    FROM    T IHE«0    TO    TWE»TFINAL 


C 


DO    600     INT«t,Ml 

CALL    RKGIINT,OI,N,W,F,LX,MX,JX) 
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PAGE    002 


XII INT>«W(l) 

Yi(  INT)«W(2) 

Kit l)«*(3> 

PI!  1)«W(<.) 
600    CONTINUE 

IF! IDUMKY.EC.l)    CO    TO    100 

IH1 1DUMMY.E0.2)    GO    TO    200 

IF(  IDUMMY.E0.3I    CO    TO    300 

IF C IDUKKY.EC.4)    CO    TO   *00 
100    Ctf.TlNUE 

DO    31    I-l,« 

XP!  1>*X1(IJ 

YP(I)«YU  II 

31  CONTINUE 
C 

C  SET    THE    INITIAL    VALUES    OF    HOMOGENEOUS    SOLUTION    til 

C 

X1(1>«0. 

YKD-1. 

R1(1>»0. 

P1I1)«0. 

IDUMMY=2 

ICHECK-2 

GC   TC   20 
200    CONTINUE 
C 
C  HOMCGENECUS    SOLUTION    (II 

c 

DO  32  I«lfHl 
XHHII«XltI) 
YHUII-Yltl) 

32  CONTINUE 

c 

C  SET     INITIAL    VALUES    OF    HOMOGENEOUS    SOLUTION    121 

C 

X1(1)«0. 

Ylt  l)«0. 

Rllll'l. 

P1I  1)»0. 

10U«MY«3 

1CHECK«2 

GO  TO  20 
300  CONTINUE 
C 

C      HOMOGENEOUS  SOLUTION  (21 
C 

DC  33  I«l,Kl 

XK2! I)=Xlt II 

YH2{1)«YIIII 

33  CONTINUE 
C 

C      INITIAL  VALUES  OF  HOMOGENEOUS  SOLUTION  (31 
C 


Xl(l>«0. 
Yl(l)«0. 
RU  1>«0. 
Pl(l)«l. 
10UMP.Y-4 
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PACE    003 


ICHECK«2 

CD    TO    20 
400    CONTINUE 
C 

C  HOMOGENEOUS    SOLUTION    C3) 

C 

DO    34    l-i.HI 

XH3U>»X1(I> 

YH3(I>«YI<11 
34    CONTINUE 
C 
C  CALL    THE    SUBROUTINE    TO    EVALUATE    THE     INTEGRATION    CONSTANTS 

C 

CALL    TUNC  I XP , YP ,XH1 , XK2 , XH3, YH1 , YH2 ,  YH  3  ,  BS  ,A ! 
C 

C  FORM   THE    GENERAL    SOLUTION 

C 

IFCA(2).LE.O.)     At2)-l. 

IF( AI2) .Gfc.10.)    A(2>*10. 

IF(A[3).LE.O.)    A<3>«1. 

IFlAI3l.GE.10.)    A{3)  =  10. 

DO   60    I«l,Kl 

XI  (I  >«XP<  I)*Atl)*XHlin*A(2)»XH2(I)*A(3>*XK3(  IJ 

Ylt  I)*Y?f  1)*A(1J*YM1(I)*A(  2)«YH2(l)t-Al3)*YH3l  I J 

60  CONTINUE 

WRITEI6.63)     NN,    At2),AC3)     ,X1  1 1 1  >  ,X1  (2  I  )  ,X  1(3  I )  ,XH  41)  ,Xl  I  5  I) , 
»XH61),XII71).XH81),XIIS1),XI{101> 
tFINK.EC.ll    GC   TO   61 
DR«RU1>-R0U) 
DP«Plll>-P0Cl> 
IF( I A3S(DR).LE.£PS) .AND.  IABS (DP) .LE.EPS) »    GO    TO    68 

61  CONTINUE 

DC    73    I»l,Ml 

xom«xui) 

Y0(  I J-Y1I  tl 

73    CONTINUE 

R0(1)«A12) 

P0I1)«A13) 

NN«NN«-l 

IFIWN-15)     90,90,68 
68    CONTINUE 

WRITE16.3I 

KA«KA«-l 

IF1KA-1)    94,94,95 
9  5    CONTINUE 

KK»KK*1 

IFIKK-1)    97,97,92 
9  2    CONTINUE 

MRITEI6.5I 

STOP 

END 
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LEVEL  21.7  (  JAN  T3  I 


OS/3  60   FORTRAN  h 


DATE   73.1  15/^3. (Jl.lS 


C 

c 
c 


FOURTH  ORDER  RUNGE-KUTTA  GILL  METHOD 

SUCRCUTINE  kKGI INT , CT , Ni Y r F ,L»H» Jl 

DIMENSION  DY(4)  ,Y(4)  ,F(28I 

T»(  1M-1)»0T 

1F1  INT.GT.ll  GO  TO  450 
410  L-3 

H«0 
450  CONTINUE 

GD  TO  (100.110. 3001, L 

100  GG  TG  (101,1101 ,IG 

101  J«l 
l«2 

DO  106  K«1,N 
rl»K  +  3»N 
K2«K1«-N 
K3*K*N 
FIKD'YCKl 
F(K3)«F(Kl) 
106  F(K.2>»DY(K) 
GO  TC  406 

110  DO  140  K«l,N 
Kl«K 

K2«K*5»N 
K3»IC2«-N 
K4»K*N 

GO  TO  (111,  112, 113. 114), J 

111  FIK1)«DY(K)*0T 
Y1K)»F(K4>*0.5*F(K11 
GO  TO  140 

112  F!K2>«0Y(K}»0T 
GO  TO  124 

113  F(K3)«0Y(K)»DT 
GC  TC  134 

114  Y(K)»F(K4)*(F(Kl)*2.«(F(K2)»FIK31)+DY(lO*DT)/t. 

GO  TC  140 
124  Y1K!«0.5»FIIC2) 

Y(K)»Y(K)«-FtK4) 

GO  TC  140 
134  Y(K)«F(K4I*F(K31 
140  CONTINUE 

GC  TC  1170.180,170,1801, J 
170  T»T*0.5»0T 
ISO  J-J+l 

IFtJ-41    404,404,299 
295    H«l 

GO    TO    406 
300    IG-1 

GO    TO    405 

404  IG-2 

405  L*l 

436    CONTINUE 

IFCP-ll    475,410,475 
475    GO    TO    (500,600,6001 ,L 
500    CALL    DFY(Y.OY) 
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GO    TC   450 
600    RETURK 
END 
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LEVEL    21.7    t    JAN    73    I  CS/360      FORTRAN   H  DATE       7S.1  15/2  5.21.50 


SUBROUTINE    DFYU.F) 

DIMENSION    21*1 ,FU> ,XO 1131) ,  YO (101) ,R0C1),P3U) 

COMMON/ ICCNTR/ I .METHOD 

CCMMCN/FORHER/X0.Y0,R0,P0 

F11W12) 

F13>»0. 

F(4)*0. 

lFlMETHOD.EQ.l)    GO    TO    10 

IF!  METHOD. EC. 2)    GO    TO    20 
C 
C  PARTICULAR    SOLUTION 

10    F{2>-P3ll)»Zl2)»2.*P3tl>»R0U)*X3(  t  J»2tl  l*P0U»»X0U>»»2*Zt3>* 

»(Y0tI)+R011)"XO(I)**2)»il*>-C3.«P0U)«k:il>»XCU>*«2+PCU)*YOll)> 

GO  TO  30 
C 
C      HOMOGENEOUS  SOLUTION 


C 


2  0  F(2)-POtl>*2l2)*2.»P0ll)*RCll>*X0(n*Ztl)«-PO(l)«X0n)««2»i(3>* 

•  1Y01  I>*R0UI*XO<I>**2>*£(*» 
30  CONTINUE 
RETURN 
ENO 


85 


LEVEL    21.7    I    JAN    73    »  CS/360       FORTRAN    H  DATE       II.1U/U.U.M 


C  SOLVE   THE    EQUATIONS    TO    UCTA1N    THE    INTUSRATICN    CCKSTASTS 

SUBROUTINE    FUKCIXP, YP.XH1, XH2.XH3 ,YH1 , Yh2, YH3 ,BS ,AA) 

DIMENSION    XPUOI).XHIUG1>,XH2110U.XH3U01     ,5S(10     .A(2;2),AA13>, 

»6i2>,YHHion,YHzeion.Yh3tioi>,rPtion,cuioj,c2UO>,Q3UO) 

SUM1-0. 

SUK2«0. 

SUM 3*0. 

SUM^»0. 

SUM5«0. 

SUN  6*0* 

DC    100    1-1.  10 

AX«XhltI»10U)/YHM10l) 

Oil  I)«XPI  !*1D*1>-AX»YP(101) 

C2(  I)«XH2  1  !»10*1>-AX»YH21101) 

03 1  1  )«XH3UMOH>-AX»YH3U01) 
100    CONTINUE 

DC    105    I«1»10 

Sim»SUiU+w2(I)*Q2tI) 

SUK2«SUM2*C3lI>*C2t I) 

SUM}*SUM3-C2U>*(01U>-BS(  II) 

SJM««SUM<.f03(  I)  •02(1) 

SU!-.5«SUM5  +  u3in*Q3Ul 

SUM6«SUMb-C3lI>*IClU>-BS( III 

10  5    CONTINUE 

A(l.ll«SUHl 

A(1,2)«SUM2 

B(i)»SUM3 

At2.1)«SUM<t 

A(2,2)«SUK5 

B(2)»SUH6 

C  SOLVE    THE    T*C    SIMULTANEOUS    EQUATIONS 

AM2>-IAl2.2>»B<l>-AU.2)»El2>>/IAU,n*A(2,2>-All.2>»AC2,l)) 

AA{3)»(6(ll-A(l.l)*AAI2)l/A(1.2>  ....... 

AA{n.-tYPU01)»AAl2)»YH2(10ll*AAl3)»YH3(lCll)/YHH10l) 

RETURN 
ENO 


PO 


1.0 
1.0 


GENEfcAL    SCLUTICNS    AT    CORRESPONDING    CATA    POINTS 
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1 

1.00000 

10.00000 

0.75  0. 

z 

1.73550 

6.83985 

0.74  0. 

I 

1.87826 

7.64583 

0.74  0. 

4 

1.97978 

5.88460 

0.74  0. 

5 

2.00302 

5.96851 

0.74  0. 

6 

1.99816 

6.02124 

0.7*  0, 

7 

1.99789 

6.02706 

0.74  0, 

6 

1 . 99647 

6.04626 

C.74  0 

9 

1.99659 

6.04449 

0.74  0 

10 

1.  998*2 

6.01924 

0.74  0 

11 

1.995*9 

6.05949 

0.74  0 

12 

1.99937 

6.00916 

0.74  0 

13 

1.99558 

6.05777 

0.74  0 

14 

1.99649 

6.04573 

0.74  3 

IS 

1.99431 

6.07467 

0.74  0 

.69  0.64  0.61  0.60  0..4I  0.64  0.70  0.76  0.90 
.67  0.60  0.55  C.51  C.47  0.44  0.41  0.40  0 .40 
.67  C.60  0.55  C.51  C.47  C.44  0.41  0.39  0.39 
.67  0.60  0.55  0.51  0.47  0>44  J .41  0.4J  0.39 
.67  C.60  0.55  0.51  C.47  C.44  0.41  0.39  0.39 
.67  0.60  0.55  0.51  0.47  0.44  0.41  0.39  0.39 
.67  C.60  0.55  0.51  0.47  C.44  0.41  0.39  0.39 
.67  C.60  C.55  0.51  0.47  0.44  0.41  0.39  0.39 
.67  0.60  0.55  0.51  0.47  0.44  0.41  0.39  0.39 
.67  0.60  0.55  C.J1  0.47  C.44  0.41  0.39  0.39 
.67  0.60  0.55  0.51  0.47  0.-.4  0.41  0.39  0.39 
.67  C.60  0.55  C.51  C.47  C.44  0.41  0.39  0.39 
.o7  0.60  0.55  0.51  0.47  0.44  0.41  0.39  0.39 
.67  0.60  0.55  C.51  C.47  0.44  0.41  0.39  0.39 
.67    C.60    0.55    C.51    C.47    C.44    0.41    C.39    0.39 
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LEVEL    21.7    {    JAN    73    1  CS/360      FORTRAN    H  DATE       73. 116/22. 08. 09 


C 

C  THE    PARAMETERS    IN    THE    NONLINEAR    DIFFUSION    EQUATION    FOR    ORYING 

C  OF    CONCRETE    ARE    DETERMINED    8Y    THE    GENERALISED    NEWTCN-RAPHSON 

C  KETHCD    (Oft    THE    QUAS ILINEAR IZATICN    METHOD).    THE    CRANK-NICOLSON 

C  TECHNICUE    HAS    BEEN    USED    FOR    NUMERICAL     INTEGRATION.    USE    OF    THE 

C  THOMAS    ALGORITHM    HAS    ALSO    BEEN    MADE. 

C 

REAL    L 

DIMENSION    Y0(  101,10)  ,YH  101,10)  ,YH1<  101,10)  ,YH2U01,10)  ,YH3U01,10 
*),YP(101,10),A10(1),A20{1)  ,A30U)  ,  Al  1 1  11  ,A21  (  1 )  ,  A3  1  (  1)  , 
*A(3),V:(10),L{10;,D(10),U<10),B[10),8S19),OP{9),QU9),Q2<9),Q319) 

CCMMCN/FDRMER/YO,A10,A20,A30,MO,M2,M3,M 

CCMMCN/ICCNTR/INT.1CHECK 

1  FORMAT (nFlO.S, 615) 

2  FORMAT {10X.9F10. h) 
4    FORMAT < 10X,3F3.4) 

6  FORMAT (515) 

7  FCR»AT(F1C.5J 
9    FDRMAT(3I5) 

13    FCRMAT13F10.6) 

63    F0?.MAT(5X,I4,6X.3F12.'t) 

RE A  0(5, 1)     DX, Tr, EPS, XE,N, Ml, MO, M2,K3,ND ATA 

READ(5,6)     1X,IY,U,  IX1.IY1 

RcAC(5,6)     2A, IS, IC, ID, IE 

DO    11     1=1,NDATA 

R£AD(5,7)     BS(I) 

11  CONTINUE 

RE6D(5,9)      INN,IX2,IY2 
READ(5,13)    DT1,0T2,DT3 
REAC(5,13)    DT1C,DT20,DT3C 
READ(5,13)     A1INT,A2INT,A31NT 
KRITF(&,4)    A1INT,42INT,A3INT 
WRITE(6,2i     (BS( i), I=1,NDATA) 
T  =  0. 

DO    12    I=1,N 
S=T+1. 

DO    120    J»1,K2 

YOJ I»J)*I l./S)**(0.CC294  2*(ALOG(S) )**2 .26) *( DX*{ J-l J )** 10. 12882- 
'  *0.03179*I  l.-DX*(J-l)  )*»2) 
YOU  ,M3)=YC(I,M0) 
120    CONTINUE 

IFt I.GE.1.AND.I.LE.IX1)    CT=DT10 
IF(  I .GE. 1X2. AND. I. L E.I  YD     DT-DT20 
IF(  I.GE.IY2.AND.I.LE.M1)    DT=DT30 
T=T+DT 

12  CONTINUE 

WRITE (6,2)    Y0( IX,IA),YO{  IX,IB),Y0(  IX, IC) ,Y01  IX,  ID)  ,Y0(  IY,IA) ,Y0{ IY 
*,1PJ,Y0(IY,IC),Y0(IY,ID),Y0(IZ,IA1 

A10(1)=A1INT 

A20( 1)=A2INT 

A30( 1)=A3INT 

NN=1 
90    CONTINUE 
C 

C  SET     INITIAL    VALUES    OF    PARTICULAR    SOLUTION 

C 
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00    121     1=2. M3 
Y1(1,I)=XE 

121  CONTINUE 
Yl( 1,1)=0. 
Allll>=0. 
A21 tl)=0. 
A3H  1>=0. 
IDUKMY=1 
1CHECK=1 

20   CONTINUE 

DC    122     !=1,M3 
WU  )=Y1  (  1 » I ) 

122  CONTINUE 
T=0. 

C  INTEGRATE    THE    DIFFERENTIAL    EQUATION    FROM    TIME=C    TO    TIME=TF1NAL    BY 

C  THF    CRANK    NICOLSON    METHOC,    TAKING    FINITE    DIFFERENCES    FOR    BOTH 

C  SPACE    AND    TIME    DERIVATIVES 

C 

DO    600    INT*1»K1 

IF(  INT. GE. LAND. INT. LE. IX)     DT*DT 1 

IF(  I  NT. GE. I XI. AND. INT.LE.IY)    0T=DT2 

IM INT.GE.IY1.AND.INT.LE-M1)     DT=DT3 

R=DT/DX**2 

CALL    HATRIXlRtK,L»C,Uf B, Al 1 . A21 , A3  1 ) 

DO    123    J=1,M3 

Yl(  INT  ,J)=WU) 

123  CONTINUE 
630  CONTINUE 


IFI IDUMMY.EC.l)  GO  TO 

100 

IF( IDUMHY.EQ.2)  GO  TO 

20  0 

IFI IDUXHY.EQ.3)  GO  TO 

300 

IF( IDUMMY.E0.4)  GO  TO 

400 

PAkTICULAR  SOLUTION 

c 

C 

c 

100  CONTINUE 

DC  31  1=1, Ml 

DC  124  J=1,M3 

Y?<  I,J)=Y1CI,J) 
124    CONTINUE 
31    CONTINUE 

QPl 1  >=YP< IX, IA) 

CP(2)=YP(IX,IB) 

0P( 3)=YP( IX, IC) 

QP(4)*YPI IX, ID) 

QP(5)=YP( IY.IA) 

QP{6)=YPt IY,1B) 

0?(7)=YP{ IY, IC) 

CP(S)=YP( IY,1D) 

QP(9)=YP( IZ.IA) 

hRITE(6,2>     (CPI  I), I  =  1,NDATA) 
C 

C  SET    INITIAL    VALUES    OF    HOMOGENEOUS    SOLUTION    (1) 

C 


DO    125    1=1, M3 
Yl(  1, I )  =  0. 

125    CCNTINUE 
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All(l)*l. 

A2K  1)=0. 

A31U)=0. 

IDUMMY=2 

ICHECK=2 

GO    TC    20 
C 

C  HOMOGENEOUS    SOLUTION    (1) 

C 

2J0    CONTINUE 

DC    32    1=1, Ml 

DO    126    J=1,M3 

YHK  I,  J)  =  YKI  ,J> 

126  CC.NT  INUE 

32  CONTINUE 

OK  1)  =  YH1<  IX,  IAJ 

CH2)«YH1(  IX,  IB) 

Q1(3)«YH1CIX(IC) 

0K4)=YH1( IX, ID) 

CK5)=YH1 { IY.IA) 

01(6)=  YHK  IY,1B) 

Cl(7)=YHl (IY.IC) 

SK8)=YH1  1IY.ID) 

QK9)  =  YHKIZ,IA) 

WRITE16.2)     (CHI),  I  =  1,NDATA) 
C 

C  SET    INITIAL    VALUES    OF    HOMOGENEOUS    SOLUTION    12) 

C 

DO    127    J=1,M3 

YK  1,J)=0. 

127  CONTINUE 
A1K  1)=0. 
A21 (i)=l. 
A3K  1)=0 
IDUMMY=3 
ICHECK=2 
GO    TC    20 

C 

C      HOMOGENEOUS  SOLUTION  (2) 

C 

300  CONTINUE 

DO  33  1=1 ,M1 

DO  128  J=1,M3 

YH2II, J)=Y1(I,J) 

128  CONTINUE 

33  CONTINUE 

C2( 1 )=YH2(IX,1A) 

C2(2)=YH2l IX, IS) 

02(3)=YH2( IX, 1C) 

Q2(4)*YH2tIX.lD) 

Q2(5)=YH2(  IY,  U) 

C2I6)=YH2CY,  IB) 

Q2(7)=YH2(IY,IC) 

C2(B)=YH2(IY, 10) 

G2(9)=YH2(IZ,  1A) 

WRITEI6.2)     (Q21I) ,I=1,NDATA) 
C 
C  INITIAL    VALUES    OF    HOMOGENEOUS    SOLUTION     (3) 
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DC  129  J=1,M3 

Y1(1,J)=0 
125  CONTINUE 

AIM  1>=0. 

«l(ll-0. 

A31!  1)=1. 

IDUMMY=4 

ICHECK=2 

GO    TO    20 
C 

C  HOMOGENEOUS    SOLUTION    13) 

C 

40  0    CONTINUE 

DO    34     1=1 ,K1 

DO    130    J=1,M3 

Yri3( I,J)=Yl(I,J) 
1.30    CONTINUE 
34    CONTINUE 

03!  1 >=YH3[IX,1A) 

03(2)  =  YH3(  IX, IB) 

C3m  =  YH3<  IX,  10) 

C3(4)=YH3l IX, ID) 

03(5)=YH3UY,IA) 

C3(fc)«YH3< IY, IB) 

Cli(  7}=YH3(IY,1C) 

Oil  S)=YH2( IY, ID) 

C3(9)=YH3l 12, IA) 

WRITEI6.2)     (03( I) ,I=1,NDATA) 
C 

C  CALL    SUBROUTINE    TO    FIND    INTEGRATION    CONSTANTS 

C 

CALL    FUNC(CP,Cl,C2,Q3,BS,A> 

KRITEI6.63)     NN,A(1 ) ,A(2) fA(3) 
C 

C  GENERAL    SOLUTIONS    FOR    Yl (T 1M E, SPACE) 

C 


DO    60    I«1,K1 

DC    131     J=2,M2 

Yll  I  ,J)=YP( 1,J)*A(1)*YH1 (I ,J)+A(2)*YH2(  I  , J > +A £3) *YH3 (I , J) 

131  CONTINUE 
60  CONTINUE 

WRI  TEC  6.21  Yl(  I  X.I  A)  ,Y1(IX,IB)  ,Y1<  IX.IC)  ,Y1(IX,ID)  .YKIY.IA)  ,Yl(IY 
*,I5),Y1{IY,IC),Y1(IY,ID),Y1(U,IA) 

DA1«   All)-AlO(l) 

DA2  =   A(2)-A20(l) 

DA3=   A(3)-A30(l) 

IF!  (  ABS  I  DAD.LE.EPS).  AND.  (ABS(DA2)  .LE.EPS)  .ANC.  J  ABS  (  CA3  )  .LE.EPS  )  ) 
*  GO  TO  66 

DO  73  1=1 ,M1 

DO  13?  J=1,M3 

Y0(  I,J)=Y1(  I, J) 

132  CUNTINUE 
73  CONTINUE 

AlOtl)=A(l) 
A20(1)=A(2) 
A301 1)=A(3) 
NN=NN*1 
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63 


CONTINUE 

STOP 

END 


90,90,66 
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LEVEL  21.7  (  JAN  73  )  CS/360   FORTRAN  H  UATE   78.116/22.09.06 

SU8ROUT I NE  MATR I  X ( K , Z , L , C , U , b , Al N , A2N ,  A3 M 

REAL  L 

DIMENSION    Zf  10)  ,Y0(101  ,10)  ,A1N(1 ),  A2NI  l),A3N(l).Li  lOltlillO)  ,D(  10  J, 

*  R(  10)  ,A1( 1) ,A2(1) ,A3U) 
CCMfCN/FORMER/YCfAl,A2»A3f«0»M2i«3i« 
COMMON/ ICONTR/1 .METHOD 

Z(1)«0. 
ZIK3)=Z(M0) 
DC    111    J=2,K2 

L(J)    =    -R*CUli 1) ,A2ll) ,A3t l),YOt I ,J) ,Y0(1*1, J)) /2. 
D(J)    =    R*{1./R*CU1 U) ,A2(  1) ,A3(  I) ,Y0( I, J) ,Y0<  1*1, J)  )♦ 
»    CY1A11  1)  ,A2(1)  ,A3( 1) ,YO(I  ,J)  ,YO(I+l  ,  J  ))/*..*(  YO  C l+i,J)-YO( I+i,J-l) 

*  +  YO(l,  J)-YO£  I.J-1)  )-CYY(Al(l),A2(  1 )  ,  A3(  I )  ,  Y0(  I  ,  J  )  ,  YO  [  H-l ,  J)  )  /8.* 

*  (YOtl  +  l,J  +  l)-YO(I*l,J)+YO(l,J  +  l  >-Y0U»J))**2) 

U(J)    =    -r.*(CUi(l)  ,A2(  1)  ,A3(i>  ,YOl  I,J>  ,Y0(  I  +  l,J)>/2.* 

*  CY(A1(  I)  ,42(1)  ,  A3(  l),YOlI,J>fY0UU,J))/4.*(YOf  I  ♦•  I ,  J*  1  1-YOU  +  l  .  J) 

*  *YO(I ,J+1)-Y01I,J) ) ) 
111 


CCNT1NUE 

L(2)=0 

LtV.2)  =  UIM2)*-L(y,2) 

U(M?)=0 

IF (METHOD. 

EC.l) 

GO 

TC 

10 

IFlMETriOD. 

EC. 2) 

GO 

TO 

20 

PARTICULAR 

SOLUTIONS 

c 
c 
c 

10    CCI\T!NUE 

DO    101    J=2,M2 

BUI    =       2(J)       +  R/2.»(Y0U+l,  J+1)-2.*Y0(  I*l,J)+YCU*l,J-lJ«-YOU  ,J*l 

*  )-2.*Y0{I,J>+YO(I,J-l>)*(CllAl<l)iA2(l),A3U>,YOll,J),YJUU,J)) 

*  *IA1N(1)-Al(l))+C2(A1(1),A2(1),A3(1),YC(I  ,J)  ,  YO  11  + 1  ,  J  )  )  *  I  A2M  1 ) -A 

*  2(1)  )+C3(  Aid  ),A2(1  ),A3(  1)  ,Y0(  I,J)  ,Y0(  1*1,  JJ  )*(A3N(  D-A31  1)  )+CY(  Al 

*  (l),A2ll),A3lll,Y0(I,J),Y3(I+l,J>  )  /2.  *  (  Z(  J  J  -Y0<  I*  I ,  J»  -YO  (  1,J)  D* 
*R/<i.*(YCl  1  +  1,  J  +  D-YOI  I*l,J)  +  YO(I  ,J+i)-Y2(I  ,J>  >**2MCY1  (AH  1)  ,A2(  1) 
*,A3(  1)  .  Y0<  I,J)  ,Y0(  IU,J)  >*(A1N(1  )-AlU  )  >+CY2(/U(l)  ,A2ll),43l  1)  , 
*Y0<  !  ,J)  |yOU*1,  J))*U2M  U-A21  I)  )+CY3(  Al  (1)  ,A2tl).  A3U)  ,YOU»J)  r 
*Y0(  J  +  l  ,  J)  )*(A3N(1)-A3<  1)  >  +    CYY(  All  1 )  ,  A2  (  1)  ,  A3(  1 )  ,  YC(  1 ,  J)  ,Y0  I  1+ 1  ,  J) 
*)/2.*(ZlJ>-Y3lH-l,J>-Y0U,J>  ))«-RM.*CY(Al(  1 )  ,  A2  I  1 )  ,  A31  1 )  .YOU  ,  J)  , 

»yo(i+i,j))»(yo(I+i,j*i>-yo{i*i,j)*yo(i,j+i)-y:ii,jj)*{Z(j*1)-z{j)) 

*  +R/2.*C< Al{l),A2( 1) ,A3( 1) ,Y0( I  *  J] ,Y0(  1*1, J)  )  «  IZ I J ♦  1  >-2.*Z < J)«-Z( J- 
*1)1 

101  CONTINUE 
GC  TC  30 
C 

C      HOMOGENEOUS  SOLUTIONS 
C 

20  CONTINUE 

DC  102  J=2,M2 

B(J)  =   ZtJ)   +  R/2.*(YOU  +  l,J+l)-2.*YJU  +  l,J)+Y3(I*l,J-l)*Y0U.J  +  l 

*  )-2.*Y0t  I,J)  +  Y0U,J-1>)*(CIU1U).A2U),A2U>.Y0U,J),Y0II«-1,J)) 

*  *(A1N(1))        <-C2(Al  (1  ),  A211).  A3U),Y0(I,J)  ,Y0(I  +  1.J)  >*(  A2N(  1)} 

*  +  C31A11  1)  ,A2(  1)  ,A3(  1)  ,YGU  ,J>  ,Y0U*1,  J))*iA3N(l))        +CYIA1 

*  (l),A2tl),A3(l),Y0(l,J),Y0{I+l,J))/2.*IZ(J)))+ 

*KM.*(YOI  ]  +  l,J+l)-YJU+l,J)+YO{l,J  +  l>-YO(I,J)  )**2*(CY11A1(  1),A2(  1) 
*,A3U)  ,YOU.J),YO{  I  +  i.J)  )*<A1NU)  )  +C Y2 i All  1)  , A2( I » ,A3 (1) , 

*YOl  I.J)  fYOCUlt  JJ)*(A2N(  1)  I  *CY3(AK1),A2(1),A3(1),Y0(1,J), 
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*YO(I*li  J)  )*(A3N(1)  )  +    CYY1AK  IJ  ,A2(  1)  ,A3(  1)  iYC(  I,  JJ  fYO(I  +  l,J) 

*}/2.*(Z(J>)>  +R/4.*CY(AlU),i2{lJ,A3(U,Y0lI,J)f 

*Y0( I+l.J) )*(YO(I+liJ+l)-VO(;+l,J)+YO(I ,J+1)-YCII,J>)*J2(J*1)-2(J)> 

*  +R/2.*CiAl(l),A2ll),A3(  1)  ,Y0(  I ,  J  )  ,  Y0(  1+ 1,  J  )  )  *  (Z  U  + 1  >-2.*Z  (J  >+Z(  J- 
*1J  ) 

102  CONTINUE 
30  CONTINUE 

CALL  TRI0AG(2,H2»LfCfU»B,Z) 

RETURN 

END 

FUNCTION  CY3{A1,A2,A3, YitY21 

CY3=  (l.-Al)*A2**A3*U  .- (0  .5*(  Y1+Y2  )  )*( A3-1 . ) + (A3- 1. )*( A3-2 . )/2.* 
*( 0.5*1  Y1+Y2) )**2-(A3-l.)*{ A3-2.)*( A3-3.)/6.*(0.5*t Y1+Y2) )**3>*( 
*(l.+A3*(AL0G{A2>-{0.5*(Yl+Y2})-  I  0. 5*( Y I  +  Y2J  ) **2  11.-    (C.5MY1  +  Y2) 
*)**3  /3.-  (0.5*(YH-Y2n**4  Ik.-    (G.5*(Yi+Y2) >**5  / 5 .-  (J.5*(Y1+Y2 
*) )**fc  /a.) )+ 

*U.-A3*(ALCG<  A2)-(0.5*1Y1  +  Y2)  )-  (  0  ,5*(  Y  1  +  Y2)  )  **2  /I.-    <C.5*IYH-Y2) 
*)**3  /3.-  ( 0.5*(Y1+Y2) )**4  M.-  i 0 .5*< Y 1 +Y 2 ) ) **5  / 5 .-  10.5*(Y1+Y2 
*) )**fc  /6.  )  )  *A2**A3* 

*[ l.-A3*(0.5*(Yl+Y2) )+A3*(A3-l. >/2 .* ( 0. 5* (Y 1+Y2 ) )**2-A3* (A3- 1. ) * 
*(A3-2. J/6.*(0.5*(Y1+Y2) J**3) )/ 

*  (l.+A2**A3*(l.-A3*(0.5*(Yl  +  Y2) )+A3*U3-l. ) II  .*( C. 5*1 Y1+Y2) )**2- 
*A3*(A3-1.)*(  A3-2.)'/6.*(0.5*(Yl+Y2)  )**3)  }**3 

RETURN 
END 


FUNCTION  CY2( Al ,A2, A3, Y1,Y2) 
CY2=  ( l.-AL)*A3**2*A2*(A3-l.)* 

*  (  l.-(  A3-1.  )*(C.5*(  Y1+Y2)  )  +(  A3-1  .  )  *(  A3-2.)  /2.-MG.5*  (Y1  +  Y2)  )**2- 

*  (A3-1 . J*(A3-2. )*( A3-3.)/6.*(0.5*{Yl+Y2J )»*3)*(1 .-A2**A3* 

*  ( l.-A3*( 0.5*(Y1  +  Y2)  )  +  A3*<  A3-I.J /2.*(0.5*( Y1+Y2) >**2-A3*U3-i. i* 

*  (A3-2. J/6.*(C.5*{Y1+Y2>)**3) >/ 

*  (  l.+A2**A3*<  l.-A3*(0.S*(Yl+Y2  >  )  *-A3*  (  A3- 1.  )/2  .* ( 0. 5* (Y 1+Y2 )  }**2- 
*A3*( A3-l.)*U3-2. )/6.*( 0.5*1 YI+Y2J )**3) )**3 

RETURN 
END 
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FUNCTION    C(Al»A2|A3iYl,Y2l 

*  ll.*A2*»A3*« l.-A3*(0.5*(Yl+Y2))*A3*lA3-l.J/2-*(0.5*(Yt*Y2))**2- 

*A3*(A3-l.)*(A3-2.)/6.*{0.5*(Yl*Y2) )**3)» 
RETURN 
END 

FUNCTION    CYUl,A2,  A3,Y1,Y2> 
CY    =    A3*(  l.-Al)*A2**A3* 

*  li.-(  A3-l.)*(0.5*t.YUY2>>-H»3-l.l*U3-2.)/2.*ia.5*CYH-Y21  )**2- 

*  (A3-1. )*(A3-2.)*( A3-3. I /6.*(3 .5*1 Yl +Y2) )**3 J / 

*  t  l.*A2**A3*J  l.-A3*(0.5*(Yl+Y2)»«-A3#(\3-l.)/2.*{0.5*CYUY2)  >**2- 
*A3*lA3-l. )*(A3-2.)/o.*(0.5*(Yl+Y2) >**31 ) **2 

RETURN 
END 

^UNCTION    CYY(A1,A2,A3,Y1,Y2) 

CYY=    A  2** A  3* A3* (  l.-Al)*(l.-(  A^-2.)*t0.5*(YlfY2))+U3-2.)*{Al-3.) 
*/2.*(0.!J*{YH-Y2)  )**2-i  A3-2.)*l  A3-3.)*(  A  3-4.  )  It.*  (0.  i*(  YI+Y2)  )**3> 
**IA2**A3*< A3+1. >*(  1 .-A3*t0.5*IYl  +Y2)  \* A2*t A3-  1 . ) II . * ( 0 . 5*1 YI+Y2) )* 
**2-A3*( A3-l.)*(A3-2.)/6.*(0.5*(Yl  +  Y2))**3)-U3-l.) )/ 

*  (l.+A2**A3*<l.-A3*«0.5*(Yl+Y2)  )  +  \3*(  A3-  1.  )  12  .  *l  C.  5*(  YUY2  )  J**2- 
*A3*(A3-J.)*(A3-2.)/6.*(0.S*(Yl+Y21  >*»3J)**3 

RETURN 
END 

FUNCTI 
CI   *    I 

*  11  .+A 
*.i3*(  A" 

RETURN 
END 


ON    CHAl  ,A2,A3,YlfY2) 

.-I./ 

2**A3*(  l.-A»*(C.5*lYUY2)  )*A3*lA3-l.  )/2, 

-I.  >*( A3-2. j/o.*(0.5*( Y1+Y2)  )**3) ) 


*l O.S*IYl«-Y2J  )' 


FUNCTI 

r  2    -    — 

*  tl.-A 

*  (A  3-2 

*  (l.+A 
*A3*(A3 

RETURN 
ENO 

FUNCTI 
C3    =   - 

*  I  (  C  .  5* 
*-  (  (  0  .  5 

*  (I. -A 

*  (4  3-2 

*  tl.+A 
*A3*( A3 

RETURN 
ENO 


ON    C2(Al,A2,A3,Yl,Y2) 

(  1.-A1  >*A3*A2**( A3- I .  > * 

3*(0.5*(Y1*Y2) )+A3*{ A3-1 

.  )/6.*(0.5*(Yl  +  Y2J  >**3)/ 

2**A3*(  l.-A3*(0.5*(Yl+Y2)  )  ♦  43* (  A  3-  I.  >  II 

-1-  )*( A3-2.)/6.*(0.5*(Yl+Y2) )**3) )  **2 


)/2.*( 0.5*1 Y1+Y2J )**2-A3*(A3-l. >* 
,*( 0.5*(Y1*Y2J )**2- 


ON    C3(AliA2,A3,Yl,Y2) 

(  i.-AU*(AU.u<A2)-<0.5*(  Y1*Y2)  )-(  {  C.  5*(  YUY2)  )**2)/2.- 

(Y1+V2)  )**3)/3.-!(0.i*(Yl  +  Y2)  )**4)/4.-i  (0.5*<Y1+Y2)  J**5)/5< 

*( Y1+Y2) )**o)/6.  )*A2**A3* 

3*[0.5*(Y1+Y2) »+A3"M A3-1. )/2.*l0.5*{ Y1+Y2) )**2-A3*( A3-1 . )* 

. )/6.*( 0.5*1 Y1+Y21  )**3)/ 

2**A3*l  1.-A3M0.  5*(Y1  +  Y2)  )  i-  A  3*  (A  3-1.  )  II  .  *(  C.  5*  t  Y  1+Y2)  >**2- 

-1.  )*<  A3-2.1/fc.*(0.i>*{  Y1+Y2)  )**3)  )**2 


FUNCTION    CY1 (Al  ,A2 , A3rYl  ,Y2» 
CY1=    -A3*A2**A3* 

*  (1._l  A3-1.  )*(0.5*(Yl+Y2))+(  A3-I.XI  A3-2.)  /2.*I0.5*(YL  +  Y2)  )**2- 

*  (A3-l.>*iA3-2.)*<A3-3.J/6.*t'J.5*m+Y2})**3i/ 

*  (  U+A2**A3*I1.-A3*(C.5*-(Y1  +  Y2I  )  K\3*(A3-1.  >/2.*«  0.  5*  (Y1+Y2  >  )  **2- 
*A3*CA3-l.)*(A3-2.)/6.*(0.5*(Yl*Y2J 1**3) )**2 

RETURN 
END 
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LEVEL  21.7  (  JAN  73  )        CS/360   FORTRAN  H       DATE   78.116/22.09.17 


C 

C  USE    THE    THOMAS    ALGORITHM    FOR    THE    TRID1AGCNAL    EAND    MATRIX 

C 

SUBROUTINE    TRIDAGt IF,L,A,B,C,D,V) 

DIMENSION    Atl ) ,6(1 ) ,C( 11 ,D(1 J.Vt II ,BETA( 101) , GAMMA (101) 
C 

C  COMPUTE    INTERMEDIATE    ARRAYS    BETA, GAMMA 

C 

BETMIF>  =  8(  IF) 

GAMMA! IF)=D(IF)/BETA( IF) 

IFP1=IF+1 

DO    1    I=IFP1,L 

3ETA(I )=RII)-A(I)*C(I-l)/BETA(I-l) 

1  GAMMAl  I)  =  (DU)-A(I)*GAMMA{  1-1)  )    /BETA(I) 
C 

C  COMPUTE    F.  1NAL    SOLUTION    VECTOR    V 

C 

V(L  )=GAKMA(L J 

LAST=l-IF 

DO    2    K=1,LAST 

I  =  L-K 

2  V(I  )=GAMMA(I)-CU)*VUn  )/B£TA(  I) 
RETURN 

END 
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LEVEL    21.7    (    JAN    73    )  CS/360       FORTRAN    h  DATE      78.116/22.09.23 


C 

C      SOLVE  THE  EQUATIONS  TO  OBTAIN  THE  INTEGRATION  CONSTANTS 

C 

SUBROUTINE    FUNC(QP,C1,C2,Q3,BS,AA) 

DIMENSION    CP(9)fQl<9},C2{9),l23(9),BSl9),A(3i3),B(3)>AAC3> 

SUM  1  =  0. 
SUM 2=0. 
SUM3=0. 
SUM4=0. 
SUM 5=0. 
SUM6=0. 
SUW,7=0. 
SUMS=0. 
SUM  9=0. 
DO     ICO     1=1,9 

sumi=sumi«-ci  m*Qi  id 

SUM2«SUN2+QH I)*02II) 

SUM3=SUM3+01(I)*Q3(I) 

S'JM4=SUM4  +  G2U  >*C2  II  ) 

S'J'-'5=SUM5  +  C2(  I)*Q3(1  ) 

SuK6=SJMo+C3(I )*C3( I) 

SUV.7=SUM7-C1(I)*(QP(I>-BS(  I)  J 

SUM3*SU«8-C21I  >*OPl  1)-BS(  I)  ) 

SUr'!9=SUM9-C3n)*(0P(I)-BS(  I)  ) 
100  CONTINUE 

Atl  ,  1)  =  SUM1 

A(l  ,2)  =  SUM2 

A(1,3)=SU"3 

A(2,U=SUM2 

A(2,2)=SUM4 

A(2,3)=SUM5 

A(3,1)=SUM3 

A(3,2)=SUM5 

A(3,3)=SUM6 

B(  1)=SUM7 

B<2)=SUM8 

B(3)'SUM5 
C 

C  SOLVE    THE    THREE    SIMULTANEOUS    EQUATIONS 

C 

AA(3)    =  U3ll>*A(2,l)-B(2>*A(l,l)>*(A(2,2)*Al3,l)-A(3,2)*A(2,l))- 

*  lfi(2)*A(3, 1)-B(3)*A[2, l))*(A{ 1 ,  2)  »A  I  2  ,  1 )- A{  2  ,  2)  *A  1 1 ,1  J ))/ 

*  ([A(l,3)*A(2,l}-At2,3>*A(l,l))*UC2,2)*A(3,l)-A<3,2)*Al2,l))T 

*  <Al2,3)*A(3,i>-A( 3,3)*    A( 2, 1) )* I A  1 1 ,2 )* A{ 2 , 1 ) -A (2 ,2 )* A (1 , 1 ) ) ) 
AA(2)    =    (B(2)*AJ3, 1>-B(3)*A{ 2, 1)-AA( 3 ) »( A( 2, 3 >*AC 3, I )-A{ 3 , 3) *A « 2  ,1 

*  )))/lA(2,2)*A(3,l)-A{3,2)*Al2,l)) 

AA(1)    =    tbll)-AA(3)*A( 1, 3)-AA{2)*A(l,2J) /A(l,l) 

RETURN 

END  , 
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ABSTRACT 

Application  of  the  quasi! inearization  technique  for  parameter  estimation 
in  nonlinear  differential  equations  is  investigated. 

Parameter  estimation  of  the  nonlinear  differential  equation  for  a  homo- 
geneous tubular  flow  chemical  reactor  with  axial  mixing  is  presented  in  Part  I 
It  is  assumed  that  the  physical  process  can  be  represented  by  a  nonlinear 
ordinary  differential  equation  of  known  form  but  containing  unknown  para- 
meters. Linearization  of  the  problem  is  carried  through  by  the  quasi- 
linearization  technique.  The  fourth  order  Runge  Kutta  Gill  method  is  used 
for  numerical  integration. 

An  algorithm  is  devised  based  on  the  least  squares  method.  It  min- 
imizes the  error  between  the  experimental  data  and  the  model  predictions. 
Numerical  data  from  Lee  [11]  has  been  treated  as  the  experimental  data. 

Extension  of  the  application  of  the  quasi  linearization  technique  to 
nonlinear  parabolic  differential  equations  is  investigated  in  Part  II. 
The  mathematical  model  of  moisture  diffusion  in  concrete  medium  is  pre- 
sented. The  Crank  Nicolson  method  is  used  for  numerical  integration. 
Experimental  data  of  Abrams  and  Orals  [3]  is  used. 

The  method  converges  quadratically  in  the  solution  of  the  problem 
in  Part  I.  Convergence  was  not  attained  in  the  solution  of  the  problem 
in  Part  II. 


